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A  bstract 


A  method  for  the  simultaneous  structural  and  control  optimization  for  torsional 
vibration  of  a  composite  plate  which  simulates  a  spacecraft  solar  array -structure  was  developed 
in  this  study.  Included  in  the  optimization  of  the  plate  are  the  location  of  the  piezoceramic 
actuators  and  sensors  that  provide  bending/torsion  actuation,  the  control  gains,  and  the 
orientation  of  the  graphite  fibers  in  the  composite  plies.  This  research  included  the  effects  of 
using  composite  tailoring  to  promote  the  coupling  of  the  twisting-bending  mode  to  enhance  the 
damping  system.  The  plate  was  modeled  by  classical  lamination  plate  theory  using  a  linear 
elastic  strain-displacement  theory.  This  theory  was  then  incorporated  into  a  performance  index 
which  contains  both  structural  and  control  parameters  that  minimizes  torsional  and  bending 
motion  at  the  tip  of  the  plate  due  to  a  torsional  force  at  the  tip.  Along  with  the  performance 
index,  there  are  inequality  constraints  on  the  amount  of  power  to  the  actuators.  This 
performance  index  and  constraints,  along  with  upper  and  lower  bounds  on  the  design  variables, 
were  incorporated  into  an  optimization  subroutine  which  then  produced  an  optimal  design  for 
controlling  both  torsional  and  bending  vibration.  This  optimal  plate  was  fabricated  and  tested 
with  a  torsional  load  to  decide  the  validity  of  the  theory  in  improving  damping.  A  comparison 
was  made  between  the  results  of  the  theory  and  the  experimental  results  of  testing  the 
optimized  plate  and  a  baseline  plate  made  up  of  a  quasi-isotropic  lay-up.  The  frequencies  of 
the  experimental  and  the  theoretical  results  were  within  15%  of  each  other.  Also  the  damping 
factor  for  the  1st  torsional  and  2nd  torsional  modes  of  vibration  increased  significantly  for  the 
optimized  plate  versus  the  baseline  plate  which  verifies  the  basic  premise  of  the  theory. 
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SIMULTANEOUS  COMPOSITE  TAILORING  AND  BENDING  CONTROL 
OPTIMIZATION  FOR  DAMPING  THE  TORSIONAL  VIBRATION  OF  A  PLATE 


I.  Introduction 


Background 

Spacecraft  have  always  been  subjected  to  some  form  of  disturbances.  Some  of  these 
disturbances  can  be  classified  as  either  random  or  persistent.  The  random  disturbances  are 
caused  by  such  things  as  micrometeorite  impacts,  increased  solar  activity,  materia]  degradation, 
etc.  The  persistent  disturbances  can  be  caused  by  such  things  as  cryo-cooling  elements, 
vibrations  from  flexible  elements,  momentum  wheel  mass  imbalance  or  bearing  noise,  etc. 
Some  of  these  disturbances  do  not  degrade  the  performance  of  the  satellite  due  to  their  small 
amplitude  and  the  robustness  of  the  control  system.  However  these  disturbances  will  degrade 
the  performance  if  they  are  large  enough  or  if  the  pointing  requirements  of  the  satellite  are 
stringent  enough.  The  twisting-bending  motion  of  a  solar  array  panel  caused  by  either  internal 
or  external  forces  can  induce  large  vibrations  on  the  main  spacecraft  structure  called  the 
satellite  bus. 

Statement  of  the  Problem 

The  flexible  vibration  of  the  solar  arrays,  that  are  connected  to  the  satellite  bus,  are  the 
cause  of  one  such  disturbance  source  as  shown  on  the  IntelSat  satellite  in  figure  1  which  is 
found  in  reference  [1]  and  on  the  IntelSat  WWW  Home  page.  As  satellites  require  more 
power,  the  size  of  the  solar 
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array  increases  tremendously.  These  solar  arrays  are  prone  to  vibrations  due  to  their  large  size 
and  light  weight.  The  solar  arrays  can  be  subjected  to  a  multitude  of  loading  conditions  such 
as  bending,  torsion,  combined  loading,  etc.  The  bending  vibration  of  a  solar  array  can  be 
controlled  by  attaching  bending  force  actuators  and  sensors  to  the  array  at  various  locations. 
The  sensors  are  able  to  detect  mechanical  motion  and  produce  an  electrical  voltage.  The 
actuators  receive  an  electrical  signal  and  convert  it  to  a  mechanical  motion.  However,  the 
torsional  vibrations  are  not  easily  controlled.  There  are  very  few  ways  to  damp  out  torsional 
vibrations  in  a  flat  panel  and  the  ones  that  are  available  are  very  costly  in  terms  of  weight  and 
power  requirements.  This  research  centers  on  damping  out  these  unwanted  torsional  vibrations 
in  the  structure  using  a  combination  of  bending  force  actuators  and  the  structural  bending- 
torsion  phenomenon  exhibited  by  unbalanced  composite  plate  layups. 

Solution 

This  dissertation  looks  at  controlling  the  vibrations  which  could  occur  in  the  solar 
array  by  modeling  it  as  a  plate  which  will  minimize  both  torsional  and  bending  vibrations. 

This  plate  is  modeled  by  classical  laminated  plate  theory  using  a  linear  elastic  strain- 
displacement  theory  for  combined  loading.  The  plate  is  then  simultaneously  optimized  with 
respect  to  both  the  structural  and  control  system  to  isolate  the  satellite  bus  from  the 
disturbance  source  by  improving  the  damping  of  the  plate.  Improving  the  damping  of  the 
plate  results  in  suppressing  its  motion  caused  by  the  external  disturbance.  If  the  motion  of  the 
plate  is  suppressed,  then  the  satellite  bus  is  isolated  from  that  external  disturbance.  The 
control  system  performance  is  improved  by  enhancing  the  coupling  of  the  twisting  and 
bending  modes  by  altering  the  ply  layup.  Since  this  plate  is  made  up  of  composite  material, 
the  stiffness  is  optimized  by  varying  the  orientation  of  the  composite  fibers  of  each  ply.  The 
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twisting-bending  coupling  due  to  an  unbalanced  ply-layup  is  designed  to  provide  the  maximum 
energy  to  those  vibratory  modes  which  the  control  system  can  damp.  Very  few  active 
damping  systems  work  well  in  damping  out  the  torsional  motion  of  spacecraft  structures.  If 
the  energy  from  the  torsion  mode  can  be  transferred  into  a  bending  mode  by  ply  coupling, 
then  the  disturbance  will  be  damped  out  by  the  active  control  system.  In  short,  the  entire 
system  consisting  of  the  stiffness  of  the  plate,  location,  and  power  requirements  of  the 
actuators  and  sensors  is  optimized  for  maximized  performance  for  torsional  loading.  A 
comparison  is  then  made  experimentally  between  the  modified  structure  versus  the  unmodified 
structure  for  improved  performance.  The  root  mean  squared  (RMS)  motion  and  damping 
factors  of  the  optimized  structure  are  also  compared  to  the  theoretical  results  of  this  research  to 
determine  if  the  modelling  of  the  plate  is  accurate. 

In  this  dissertation,  the  historical  background  will  first  be  discussed  along  with  the 
formation  of  the  theory.  Once  the  theory  is  complete,  the  optimization  is  performed  and  from 
this  output  an  optimized  plate  is  fabricated  and  tested.  Finally  comparisons  are  made  between 
the  theoretical  predictions  and  the  experimental  results. 
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II.  Historical  Development 


Introduction 

In  this  chapter,  the  work  accomplished  by  previous  researchers,  is  discussed  as  it  . 
directly  relates  to  this  dissertation.  Also,  the  development  of  the  fundamental  composite 
equations  of  the  theory  are  discussed.  The  historic  development  is  divided  into  two  different 
sections.  The  composite  material  theoretical  equations  are  provided  in  the  first  section. 

Instead  of  dealing  with  the  historic  development  of  composite  material,  the  elasticity  equations 
are  formulated  which  will  be  used  later  in  the  theoretical  development.  These  elasticity 
equations  are  derived  from  Hooke's  Law  for  three  dimensional  stress-strain  relationships. 

The  second  section  provides  a  background  look  at  the  theory  of  composite  tailoring 
and  coupling.  This  section  will  present  historical  work  in  this  area  about  aeroelasticity  and 
current  work  which  deals  with  the  problem  defined  in  this  research.  Aeroelasticity 
predominately  deals  with  aircraft,  but  the  same  principles  can  also  be  applied  to  parts  of  the 
spacecraft  such  as  the  solar  arrays.  The  solar  arrays  sometimes  act  as  large  wings  on  a 
spacecraft  when  in  the  presence  of  the  low  Earth  orbit  environment.  Through  the  use  of 
aeroelastic  tailoring,  better  control  can  be  accomplished  of  the  solar  array  structures. 

Composite  Material 

A  composite  material  consists  of  a  reinforcement  to  provide  strength  and  stiffness  and 
a  matrix  material  that  acts  as  a  glue  to  hold  the  structure  together.  Early  examples  of 
composite  material  include  the  clay/straw  mixture  used  to  make  adobe  structures  to  steel 
reinforced  concrete  used  to  support  very  large  structures  such  as  bridges,  dams,  etc.  When  one 
speaks  of  advanced  composites,  one  usually  means  composites  that  use  stiffening  agents  such 
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as  glass,  graphite,  carbon,  boron,  and  metal  fibers;  and  matrix  materials  such  as  epoxy,  metal, 
and  carbon  versus  the  previously  mentioned  materials  and  glass  composites.  Composite 
materials  usually  have  very  good  stiffness  and  strength  characteristics  when  compared  to 
conventional  materials  especially  when  weight  is  a  factor.  The  specific  stiffness,  or  stiffness  to 
weight  ratio,  and  the  specific  strength,  or  strength  to  weight  ratio,  of  composites  are  vastly 
superior  to  conventional  materials  like  aluminum,  steel  and  even  titanium. 

All  composite  materials  whether  the  are  hand  layed-up,  filament  wound,  woven,  or 
pultruded  are  made  up  by  stacking  individual  layers,  called  laminae,  together  to  form  the 
structure.  Each  lamina  consists  of  fibers  oriented  in  the  same  direction  held  together  by  the 
matrix  material.  (See  figure  2.)  The  laminae  are  then  stacked  together,  allowing  the  fiber 
orientation  to  vary  from  one  to  another  to  form  the  laminate  or  structure.  One  of  the 
principal  advantages  of  composite  materials  is  that  each  individual  lamina  can  have  varying 
fiber  orientations  with  respect  to  one  another.  By  stacking  laminae  with  varying  fiber 
orientations,  the  properties  of  the  structure  will  change.  For  example,  to  obtain  maximum 
longitudinal  and  bending  stiffness,  all  of  the  fibers  in  each  lamina  need  to  be  oriented  in  the 
axial  direction.  The  fiber  orientation  will  decide  the  strength,  stiffness,  and  almost  every  other 
property  of  the  structure.  Figure  3  shows  how  each  ply  can  be  different  from  the  preceding 
ply  by  changing  the  orientation  of  the  fibers.  This  allows  the  structure  made  out  of  composite 
material  to  be  tailored  to  meet  any  designers  requirements.  The  plot  in  figure  4,  made  by  an 
engineer  at  the  Phillips  Laboratory  using  in-house  composite  test  data,  shows  the  design 
versatility  of  using  composite  material  with  respect  to  aluminum.  The  properties  of  aluminum 
are  set  to  unity  and  the  relative  composite  properties  are  shown.  At  least  in  the  aerospace 
industry,  the  trend  is  to  fabricate  almost  everything  out  of  composite  materials  due  to  their 
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tailorability  and  superior  properties. 

Several  assumptions  are  made  to  analyze  composite  laminae:  [2] 

(a)  The  fibers  are  homogeneous,  linearly  elastic  and  isotropic.  They  are  equally 
spaced  and  perfectly  aligned. 

(b)  The  matrix  is  homogeneous,  linearly  elastic  and  isotropic.  Hooke's  Law 
governs  the  mechanical  behavior  of  a  composite  lamina. 

Jones  [3]  and  Murray  [4]  provide  a  very  good  description  of  the  macromechanical  behavior  of 
a  lamina  and  a  laminate.  The  generalized  Hooke's  Law  for  an  anisotropic  material  can  be 
written  in  matrix  form  as: 
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where  Gj  are  the  stress  components,  [C]  is  the  stiffness  matrix,  and  £j  are  the  strain 
components.  The  stiffness  matrix  can  be  related  to  standard  engineering  constants  such  as 
Young's  Modulus  E,  Possion's  Ratio  v,  and  the  shear  modulus  G  for  each  direction.  In 
anisotropic  materials  there  are  no  planes  of  symmetry,  since  none  of  the  axes  of  the  lamina 
line  up  with  the  principal  axes  of  material  symmetry. 

In  most  composites,  two  axes  of  the  lamina  line  up  with  the  principal  axes  of  the 
material  which  form  two  planes  of  symmetry.  The  1  direction  consists  of  fiber  dominated 
properties  while  the  2  direction  consists  of  matrix  dominated  properties.  This  material  is 
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defined  to  be  an  orthotropic  material.  The  stress-strain  relationship  for  this  material  is: 
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If  a  laminate  is  thin  enough,  a  state  of  plane  stress  exists  and  is  characterized  by: 

°3  =  T23  =  *31  =  0  (3) 


This  is  shown  in  figure  5  which  also  shows  the  principle  directions  of  the  structure.  The 
stress-strain  relationship  expressed  in  equation  2  is  now  reduced  to: 
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where: 
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fts  -  cM  -  an  m 

When  the  principle  material  directions  move  to  a  different  coordinate  system  and  differ  from 
the  principle  directions  of  the  structure  as  in  the  case  of  angle  ply  layups,  the  stresses  and 
strains  need  to  undergo  a  coordinate  transformation.  (See  figure  6).  The  angle  9  is  the  angle 
from  the  x-axis  to  the  1-axis.  The  transformation  from  the  1-2-3  coordinate  system  to  the  x-y- 
z  system  is: 


where  [T]'1  is  [5]: 


cos20 


m_1= 


sin20 


sin0cos0 


sin20  -2sin0cos0 

cos20  2sin0cos0 

-sin0cos0  cos20  -  sin20 


(10) 


Using  the  coordinate  transformations,  the  relationship  between  stress  and  strain,  in  any 
arbitrary  direction,  can  be  found.  This  relationship  is: 


ai) 


where 


(12) 


The  stress-strain  relationships  for  transverse  isotropy  in  xyz  coordinates  are  then  represented  in 
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matrix  form  as: 


(13) 


in  which: 

(?u=C11cos40+2(C12+2C66)sm20cos20+C22sm40  (14) 

Qn =(Cj  t  +C22  -4C66)sin20cos20 + C12(sin40  +cos40)  (IS) 

<316=(Cu-C12-2C66)sin0cos30+(C12-C22+2C66)sin30cos0  (l6) 

Q21  =Cnsin40  +2(C12+2C66)sin20cos20  +C22cos40  (1 7) 

(p26 =(CU  -C12  -2C66)sm30cos0  +(C12 -C22  +2C66)sin0c°s3O  (18) 

Q6 6<cu+c22  ~2Cu  -2C66)sm20cos20  +C66(sin40  +cos40)  ( 19) 


Now  that  the  stress-strain  relationships  for  any  arbitrary  ply  angle  are  developed,  the 
equations  relating  strain  to  displacement  for  an  anisotropic  plate  need  to  be  derived.  These 
equations  will  be  derived  in  chapter  III  using  the  strain  relationships  derived  in  a  rectangular 
cartesian  coordinate  system.  A  linear  elastic  theory  will  be  used  to  describe  this  relationship. 
The  coordinate  system  used  in  this  development  is  shown  in  figure  7.  The  x  and  y  axis  are 
located  at  the  mid-plane  of  the  laminate.  The  displacements  U0,  V0>  and  W0  are  the  laminate 
mid-surface  displacements  in  the  x,  y,  and  z  directions  respectively. 

Composite  Tailoring 

Composite  tailoring  has  been  used  in  aircraft  and  helicopter  applications  in  the  recent 
past.  The  tailoring  of  the  composite  material  properties  such  as  the  strength,  stiffness. 
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thermal  coefficients,  and  coupling  terms  to  meet  an  objective  is  what  makes  them  very 
attractive  to  use  in  all  applications.  This  is  especially  true  in  the  area  of  aeroelastic  tailoring 
which  is  "the  incorporation  of  directional  stiffness  into  an  aircraft  structural  design  to  control 
aeroelastic  deformation,  static  or  dynamic,  in  such  a  fashion  as  to  affect  the  aerodynamic  and 
structural  performance  of  that  aircraft  in  a  beneficial  way."  (See  [6], [7], [8], [9].)  Aeroelasticity 
predominately  deals  with  aircraft,  but  the  same  principles  can  also  be  applied  to  parts  of  the 
spacecraft  such  as  the  solar  arrays.  The  solar  arrays  sometimes  act  as  large  wings  on  a 
spacecraft  when  in  the  presence  of  the  low  Earth  orbit  environment.  Through  the  use  of 
aeroelasticity,  better  control  can  be  accomplished  of  the  solar  array  structures.  In  recent  years, 
most  of  the  effort  in  aeroelastic  tailoring  has  been  conducted  in  identifying  the  effects  of  an 
anisotropic  design  on  the  deformation  coupling  of  fixed  and  swept  wing  aircraft.  (See  [6]  [7] 
[10]  [11]  [12]  [14].)  More  recently,  the  shift  in  aeroelastic  tailoring  has  turned  to  the  effects 
on  rotorcraft  and  rotorsystems.  (See  [2]  [13]  [15].)  Composite  tailoring  through  anisotropic 
design  is  ideal  for  this  application  since  helicopter  blades  operate  in  an  environment  consisting 
of  inertial,  aerodynamic  and  elastic  loadings.  The  benefits  of  aeroelastic  tailoring  for 
helicopter  blades  with  an  off-axis  layup  are  that  a  twist  will  be  induced  in  the  blade  when  a 
non-torsion  force  is  exerted  which  will  improve  the  performance  of  the  lifting  surface. 

In  this  research,  the  idea  of  composite  tailoring  will  be  extended  to  the  area  of 
controlling  large  space  structures  and  modified  slightly.  The  layup  will  remain  symmetrical  to 
prevent  warping,  but  it  will  be  unbalanced  which  will  cause  the  required  coupling  in  motion. 
The  plates  in  a  large  space  structure,  with  an  anisotropic  composite  layup,  will  be  optimized 
with  respect  to  their  directional  stiffness  and  coupling  terms  to  enhance  the  effect  of  the 
control  system  on  the  structure.  The  orientation  of  the  composite  fibers,  which  decide  these 
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coupling  and  stiffness  terms,  will  be  optimized  with  respect  to  the  active  damping.  The  axial- 
bending  and  torsion-bending  coupling  due  to  an  anisotropic,  unbalanced  ply-layup  will  be 
analyzed  to  provide  the  maximum  energy  to  those  vibratory  modes  that  the  control  system  can 
damp.  Very  few  active  damping  systems  work  well  in  damping  axial  and  torsion  motion.  If 
the  energy  from  the  axial  and  torsion  modes  can  be  transferred  into  a  bending  mode  through 
ply  coupling,  then  the  disturbance  can  be  damped  out  by  the  active  control  system. 

The  idea  of  composite  tailoring  for  structural  elements  in  a  spacecraft  has  been 
presented  only  recently  in  the  literature  by  Hwang  and  Gibson  [16],  Barrett  [17]  [18],  Belknap 
and  Kosmatka  [19],  Bronowicki  and  Diaz  [20],  Olcott  [21]  [22]  [23],  Hwang,  Hwang,  and 
Chul  [24]  [25]  [26],  and  others  [27]  [28]  [29]  [31]  [30],  Hwang  and  Gibson  [16]  presented 
the  influence  of  the  vibration  coupling  effects,  such  as  bending-twist  and  bending-extension, 
on  the  damping  of  laminated  composites.  The  only  damping  was  the  inherent  material 
damping  of  the  composite,  modelled  by  the  strain  energy  damping.  They  concluded  that  the 
coupling  effects  of  the  composite  are  dependent  upon  the  fiber  orientation  and  laminate 
geometry.  The  coupling  effects  were  found  maximized  at  fiber  orientations  of  30  degrees. 
Also,  the  coupling  effects  tended  to  increase  the  vibrational  damping  in  the  first  two  flexible 
modes. 

Barrett  [17]  [18]  discussed  the  improvement  of  damping  properties  of  axially  loaded 
structural  members  made  up  of  anisotropic  composite  laminated  and  damping  material.  He 
formulated  a  structural  theory  that  characterizes  the  behavior  of  an  axially  loaded  composite 
cylinder  and  applied  it  to  a  three  layered  structure.  This  structure  consisted  of  concentric 
layers  of  both  stiffness  and  damping  materials.  In  his  theory,  the  stiffened  cylinders  were 
treated  with  thin  wall  theory  that  implies  that  the  non-membrane  stresses  were  negligible. 
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Different  fiber  orientation  and  layups  were  studied  for  controlling  the  resonant  response  of  the 
structures.  These  design  study  comparisons  were  not  made  on  an  equal  mass  basis  which 
accounts  for  some  decrease  in  resonant  frequency  when  the  damping  material  is  added  to  the 
structures.  The  design  studies  proved  the  effectiveness  of  the  coupling  in  suppressing  and 
controlling  resonant  vibration.  Barrett  discussed  some  difficulties  associated  with  the  addition 
of  anisotropic  effects.  Anisotropy  will  result  in  shear  stresses  and  circumferential 
displacements  in  the  axial  composite  cylinder  construction.  These  shear  stresses  may  result  in 
new  modes  of  failure  such  as  delaminations  in  the  cylindrical  structure.  The  circumferential 
displacements  at  the  ends  of  the  members  will  cause  difficulties  because  they  are  usually 
constrained  in  most  boundary  conditions.  Barrett  discusses,  but  does  not  analyze,  some 
possible  solutions  to  these  difficulties.  One  solution  is  to  allow  the  stiffness  to  vary  along  the 
length  of  the  structure.  This  can  be  accomplished  by  allowing  the  thickness  to  vary  or  the  ply 
angle  to  change  along  the  length.  Barrett  concludes  that  the  anisotropic  composite 
construction  is  resistant  to  all  modes  of  vibration.  The  axially,  torsional,  and  flexural  modes 
of  vibration  are  resisted  by  the  coupling  introduced  by  an  anisotropic  composite  layup. 

Belknap  [19]  discussed  the  design  and  fabrication  of  a  thin  walled  composite  structure 
with  an  inner  and  outer  layer  of  graphite  epoxy  shells  and  a  center  layer  of  viscoelastic 
damping  material.  The  graphite  layers  were  fabricated  so  that  when  they  are  subjected  to 
either  a  bending  or  axial  load,  they  counter-rotate  to  produce  a  large  shear  area  at  the 
viscoelastic  layer.  Belknap  used  the  stiffness  characteristics  of  a  laminated  tube  to  choose  ply 
orientations  that  maximize  damping  and  structural  stiffness.  In  his  theoretical  development  he 
assumed  a  two-dimensional  strain  state  when  determining  the  displacement  functions.  A 
separation  of  variables  technique  combined  with  the  Ritz  method  was  used  to  solve  for  the 
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local  cross-section  warping  functions.  The  extensional-bending-twisting  coupling  can  be 
represented  by  applying  solely  an  axial  force.  By  varying  the  ply  angle,  Belknap  demonstrated 
that  the  coupling  is  maximized  when  the  angle  is  around  30  degrees.  He  also  presented  a  very 
detailed  fabrication  plan  and  problems  associated  with  incorporating  the  viscoelastic  material 
into  the  tube.  The  main  problem  that  he  encountered  was  in  testing  the  tube  in  the  axial 
mode.  He  said  that  special  end  fittings  need  to  be  fabricated  to  couple  the  stiffness  of  the  two 
composite  shells  with  the  viscoelastic  core  while  allowing  for  the  counter  rotation  of  the  two 
shells  at  the  tube's  end.  However,  by  conducting  a  simple  impact-hammer  modal  test  with  the 
tube  in  a  near  free-free  suspension,  the  damping  in  the  first  mode  was  increased  over  seven 
times  by  optimizing  the  ply  layers  and  incorporating  the  damping  layer. 

Bronowicki  and  Diaz  [20]  present  a  semi-analytical  finite  element  approach  using  a 
pair  of  concentric  cylindrical  shells  surrounding  a  viscoelastic  medium.  They  present  this 
theory  in  an  attempt  to  damp  extensional  motions  in  a  truss  member  commonly  found  in  a 
space  structure.  The  constitutive  equations  for  the  cylindrical  shells  are  derived  based  upon  a 
plane  stress,  linear  strain-displacement  relationship  for  an  orthotropic  material.  The 
viscoelastic  material  is,  as  stated  in  this  report,  "assumed  to  be  at  a  single  operating  frequency, 
which  allows  the  elastic  and  shear  modulus,  Ev  and  Gv,  to  be  defined  as  complex  numbers 
representing  magnitude  and  phase  of  stiffness."  The  analysis  that  is  performed  is  for  that 
single  frequency  and  can  be  performed  for  different  operating  frequencies  and  temperatures. 
The  stress-strain  relationship  is  derived  using  Hooke's  Law.  The  displacements  of  the 
viscoelastic  material  are  assumed  to  vary  linearly  throughout  the  thickness  between  the  inner 
and  outer  shell.  The  outer  and  inner  orthotropic  shells  are  coupled  to  the  viscoelastic  layer 
through  surface  tractions.  From  these  combined  equations,  only  allowing  radial  forces,  a  finite 
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element  is  defined  for  the  composite  structure.  These  finite  elements  are  combined  to  form  a 
coarse  grid  along  the  length  of  the  tube.  They  say  that  due  to  the  Bemoulli-Euler  hypothesis, 
the  axial  results  can  be  applied  to  bending. 

In  his  dissertation  and  other  papers,  Olcot  [21]  [22]  [23]  discusses  a  new  damping 
concept  which  uses  the  stress  coupling  effects  inherent  in  composite  materials  to  induce  large 
shear  strains  in  co-cured  damping  layers.  The  ply  angle  is  varied  along  the  length  (as  shown 
in  figure  8)  of  the  component  to  create  multiple  regions  of  high  shear  and  to  induce  high 
damping  loss  factors  within  the  component.  This  concept  was  analyzed  by  a  computer 
program  which  was  used  to  conduct  parametric  studies  on  the  analytical  models  of  flat 
membranes  and  cylindrical  damped  components.  Finally,  several  damped  cylinders  and  I- 
beams  were  built  to  test  the  validity  of  the  theory  and  computer  program.  The  model 
developed  by  the  theory  predicted  the  natural  vibration  frequency  and  damping  loss  factors  of 
the  cylinders  to  within  5%.  The  cylinders  built  and  tested  with  this  new  concept  have  tested 
loss  factors  of  up  to  8.5%. 

The  work  presented  in  the  papers  by  Hwang,  Hwang,  and  Chul  [25]  [26]  is  very 
similar  to  the  work  presented  in  this  dissertation.  Their  theoretical  model  and  approach  are 
presented  next.  However,  the  conclusions  of  their  paper  will  not  be  discussed  in  this  chapter. 
A  comparison  of  their  theoretical  and  experimental  work  to  the  current  approach  is  presented 
in  chapter  VI.  The  conclusions  of  the  paper  by  Hwang,  Hwang,  and  Chul  will  be  used  as  an 
independent  confirmation  of  this  dissertation  and  compared  with  the  results.  Their  basic 
premise  is  to  use  an  intelligent  system  concept  design  for  vibration  control  of  a  laminated 
composite  plate  with  piezoelectric  sensors  and  actuators  attached.  An  intelligent  system 
concept  design  is  another  name  for  smart  structures  which  is  the  addition  of  sensors,  actuators, 
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and  control  electronics  to  a  structure  to  control  the  vibrations.  Also  included  is  an  analysis  of 
altering  the  vibration  modes  of  the  composite  structure  by  tailoring  the  stiffness  and  taking 
into  account  the  bending-torsional  coupling.  They  determine  the  optimal  size  and  location  of 
the  sensor  and  actuator  and  the  optimal  ply  layer  angles  which  maximize  the  work  done  by  the 
controller.  If  the  work  done  by  the  controller  is  maximized  by  optimizing  the  structural  and 
control  parameters,  then  the  control  energy  utilized  is  reduced.  One  of  the  major  differences 
between  the  referenced  paper  and  this  dissertation  is  how  the  equations  of  motion  are 
determined.  Hwang  et.  al.  derive  the  equations  of  motion  by  using  Hamilton's  Principle  and 
then  discretize  them  by  using  the  finite  element  method.  The  kinematics  are  derived  by 
expressing  the  displacement  and  curvature  in  terms  of  nodal  displacements  using  the  shape 
functions  of  a  4  node,  12  degree  of  freedom  quadrilateral  plate  bending  element.  In  this 
dissertation  the  equations  of  motion  are  formulated  similarly  using  Hamilton's  Principle. 
However,  the  discrete  equations  of  motion  are  found  by  using  the  Rayleigh-Ritz  or  Assumed 
Modes  method.  Since  finite  element  theory  is  a  derivative  of  this  energy  methods,  the 
theoretical  results  should  compare  very  well. 

Also  in  the  paper  by  Hwang  et.al.,  from  finite  element  theory,  the  mass  matrix  M  and 
stiffness  matrix  K  are  generated  and  used  to  solve  the  free  vibration  analysis  for  the  natural 
frequencies  and  mode  shape  vector.  The  equations  of  motion  are  transformed  into  state  space 
form  through  the  use  of  modal  state  variables.  The  control  equations  are  adjoined  to  the  state 
space  equation  of  motion  by  using  negative  velocity  feedback.  The  combined  control 
equations  are  written  in  terms  of  both  the  structural  and  control  variables.  An  objective 
function  is  formed  based  upon  these  variables  whose  solution  will  maximize  the  work  done  by 
the  feedback  control  forces,  for  a  fixed  energy  input,  by  calculating  the  optimal  solution.  The 


21 


optimal  solution  is  dependent  upon  the  design  variables  which  are  the  ply  layer  angles  and 
number,  size,  location,  and  gains  of  each  sensor  and  actuator  pair.  When  the  work  done  by 
the  feedback  control  forces  is  maximized,  the  control  system  is  more  effective.  In  this  paper 
[25],  the  numerical  optimization  is  performed  using  the  method  of  feasible  directions. 
Numerical  calculations  are  performed  on  an  eight  ply  laminated  composite  plate  with  two  of 
the  plies  variable.  The  stacking  sequence  is  defined  as  [91/02/Oo/9O°]s  where  8,  and  9,  are  the 
ply  angle  design  variables.  The  optimization  is  first  conducted  for  an  orthotropic  plate  to 
determine  the  optimal  size  and  location  of  the  piezoelectric  sensors  and  actuators.  Initially,  the 
size  of  the  sensors  and  actuators  are  unconstrained  while  the  locations  are  not  allowed  to  vary. 
In  the  second  case  the  sizes  are  constrained  while  the  optimal  locations  are  obtained.  The 
optimization  is  next  conducted  to  determine  the  optimal  ply  angles  of  the  laminated  plate 
which  will  maximize  the  control  performance.  The  optimization  is  conducted  with  and 
without  the  dynamic  stability  constraints.  The  optimizations  are  not  performed  simultaneously 
or  integrated  as  they  are  done  in  this  dissertation.  Instead  they  are  performed  sequentially 
which  may  or  may  not  produce  the  optimal  result.  Sequential  optimization  may  be  the  easier 
method  because  it  reduces  the  number  of  parameters  which  are  being  searched  for  during  the 
optimization,  but  may  not  produce  the  best  result  as  shown  by  Dracopoulos  and  Oz  [32], 

They  show  that  using  simultaneous  or  integrated  optimization  produced  a  solution  which 
lowered  the  control  cost  by  294.4%  over  a  nonintegrated  design  procedure.  The  final  result  of 
the  papers  by  Hwang  [25],  [26]  and  in  this  dissertation,  is  that  the  control  system  performance 
increases,  as  measured  by  the  increased  damping  factors,  as  the  bending-torsion  coupling 
increases.  This  research  described  in  the  referenced  paper  presents  only  a  theoretical  analysis, 
but  in  this  dissertation  an  experiment  will  be  conducted  to  prove  the  feasibility  of  this 
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proposed  theory.  However,  the  results  of  Hwang  et.al.  are  very  close  to  the  results  found  in 
this  research  even  though  the  methodology  is  completely  different:  finite  element  method 
versus  Rayleigh-Ritz  approximation  and  sequential  versus  simultaneous  optimization. 
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III.  Theory 


Introduction 

This  chapter  is  separated  into  five  main  sections:  derivation  of  the  equations  of 
motion  for  the  composite  plate,  active  control  derivation,  derivation  of  the  complete  equations 
of  motion,  control  theory  and  structural/control  optimization  formulation  and  solutions  by  an 
optimization  subroutine.  In  each  section,  the  necessary  equations  are  derived  and  main 
assumptions  are  stated  for  each  particular  subject.  How  these  particular  sections  fit  together  is 
shown  in  figure  9  which  is  a  flow  chart  outlining  the  theory.  First  the  governing  composite 
equations  are  formulated  for  a  particular  strain  field.  These  equations  are  converted  from 
partial  differential  equations  to  ordinary  differential  equations  by  incorporating  assumed 
solutions  to  the  equations  of  motion.  Next,  the  active  control  system  is  augmented  to  the 
differential  equations  through  the  use  of  state  space  modelling.  The  equations  are  first  order 
differential  equations.  The  optimization  performance  index  and  constraint  equations  are 
determined  next.  These  equations  are  coded  into  a  FORTRAN  program  and  an  optimum 
solution  is  found  using  an  optimization  solver.  Numerous  cases  are  run  and  the  lowest  local 
optimum  solution  is  determined.  Based  upon  the  results  of  the  optimization,  an  experimental 
plate  is  fabricated  and  tested.  Finally  the  results  are  documented  comparing  the  theoretical 
predictions  with  the  experimental  results. 

Composite  Plate  Theoiy 

This  plate  is  (see  figure  10  for  dimensions)  modeled  by  classical  laminated  plate 
theory  using  a  linear  elastic  strain-displacement  theory  when  subjected  to  combined  loading. 
The  thickness  is  approximately  100  times  less  than  the  next  smallest  dimension,  the  width. 
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OUTLINE  OF  RESEARCH 


Figure  9  Outline  of  Research 


DIMENSIONS  OF  THE  PLATE 


Figure  10  Dimensions  of  Composite  Plate 
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The  deflections  of  the  plate  are  assumed  to  be  small.  Also  the  plate  is  assumed  loaded  only  in 
torsion  about  the  x-axis.  The  z-displacement,  w0,  is  a  function  of  both  x  and  y.  The  force- 
displacement  equations  and  equations  of  motion  will  now  be  stated  for  the  plate. 

Moment  Displacement  Equations.  The  generalized  displacement  field  (see  [3],  [33]) 
is: 


dwQ 

u(x,yj,t)  ~-z 

OX 

(20) 

< 

II 

«  J* 

(21) 

w(x,yj,t)  =  w0 

(22) 

The  strains  considered  in  this  analysis  are  derived  from  linear  elastic  theory.  The  linearized 
strain-displacement  relationships  are: 


= 

(23) 

dx 

dv 

(24) 

dy 

dv 

—  + 

du 

(25) 

dx 

dy 

Substituting  the  above  displacement  relationships  into  the  strain  equations,  the  strain- 
displacement  equations  become: 


= 

(26) 

s  = 

(27) 

= 

~2zwo  >*y 

(28) 
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Using  classical  laminated  plate  theory  [3],  the  moment-strain  relationship  is  presented  using 
the  above  strain-displacement  equations  which  will  be  used  to  assist  with  the  derivation  of  the 
equations  of  motion  is: 


-Du 

dI2 

D22 

d26 

D16 

D26 

(29) 


where  Qy  are  given  by  eqs  (14-19). 

The  above  relationships  are  only  valid  if  the  ply  layup  is  assumed  to  be  symmetric. 

Equations  of  Motion.  The  strain/displacement  relationship  developed  previously  will 
now  be  used  to  derive  the  equations  of  motion  and  boundary  conditions  for  the  composite 
plate.  The  equations  of  motion  are  derived  using  Hamilton's  Principle  [35]  -  [37]  which  is: 

h 

f6(T  -  U  -  V)dt= 0  (31) 

h 


where  T  is  the  kinetic  energy,  U  is  the  potential  or  strain  energy,  and  V  symbolizes  the  work 
done  by  conservative  or  nonconservative  external  forces.  The  symbol  5  represents  the  first 
variation  with  respect  to  the  dependent  variables. 

The  kinetic  energy  is  defined  as: 


27 


h 

6/2  a  2 


T  =  f  f  f  |p[(M2)  +  (v2)  +  (w2)]£feft<iy 


-6/2  0  _  ft  ‘ 
2 


where  p  is  the  mass  density  per  unit  volume  for  the  composite  plate  and  the  displacement 
derivatives  are  with  respect  to  time.  The  density  is  held  constant  throughout  the  laminate. 
The  displacement  derivatives  can  now  be  calculated  as: 

(m2)  =  z2w0,2x  ( 

(v2)  =  z2w0>y  ( 

(w2)  =  w2  ( 


Before  deriving  the  equation  for  the  kinetic  energy,  the  following  simplifications  are 
introduced: 


h  =  /  pdz 


/3  =  f  p z2dz 

2 


The  expression  for  the  kinetic  energy  is: 


oj  a 

T  =  i  f  f  +  hK2y  +  h^lxdxdy 


The  first  variation  of  the  kinetic  energy  can  now  be  found  by  taking  the  partial  derivatives  of 
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T  with  respect  to  the  dependent  variables  which  are  w0,  w0x,  w0  .  The  first  variation  is  now 


derived  as: 


8T  =  f  fWo’x^O’x  +  h*vyb*vy  + 


The  variation  of  the  kinetic  energy  is  then  integrated  by  parts  (see  appendix  A)  with  respect  to 
time  to  remove  the  partial  derivatives  from  the  variations.  For  simplification,  this  expression 
is  broken  out  into  four  parts:  contribution  to  the  equation  of  motion,  contribution  to  the  initial 
conditions,  contribution  to  the  x-axis  boundary  conditions,  and  contribution  to  the  y-axis 
boundary  conditions.  These  expressions  are  represented  as: 


I  bTdt  =Tl  +  t2  +  r3  +  r4 


Contribution  to  the  Temporal  Boundary  Conditions  can  be  stated  as: 


■  ‘2a  ■  ‘2  6/2 

Tl  =  /  hK;6wo\\dy  +  |/3>V6wol  I  <2* 

J  t  n  J  y  t  -un 


fi°  o 


/  /[(/A  -  hKxx  -  73*W6>Vo]  \dxdy 


Contribution  to  the  Equation  of  Motion  can  be  stated  as: 


h  6/2  a 


r2  =  f  f  f  [(J3%’xx  +  hKyy  "  h^dw^dxdydt 

fj-6/2  0 


Contribution  to  the  X-Axis  Boundary  Condition  can  be  stated  as: 
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Contribution  to  the  Y-Axis  Boundary  Condition  can  be  stated  as: 


h  a 


ffhKy  6w0 

A  0 


b/2 

|  dxdt 

-b/2 


The  strain  energy  is  defined  as: 

h 

b/2  a  2 

U  =  /  // +  °yey  +  V Jdzdxdy 
-6/2  0  _  h 
2 


(44) 


(45) 


or  after  taking  the  first  variation,  it  can  be  written  as: 

h 

bp.  a  2 

St/  =  jf  f  f  (0x8ex  +  ay8ey  +  ^dzdxdy  (46) 

-bp  O  h 
2 


This  expression  can  be  rewritten  after  substituting  the  terms  for  the  displacements  and 
integrating  the  z  dependence  out  of  the  expressions.  This  equation  for  the  potential  energy  is: 


6/2  a 


5  U  =  f  O’xx  +  MybwWyy  +  2Mxy8wQ,xy]dxdy 

-6/20 


(47) 


The  variation  of  the  potential  energy  is  then  integrated  by  parts  twice  with  respect  to  x  and  y 
to  remove  the  partial  derivatives  on  the  variations.  Again  for  simplification,  this  expression  is 


broken  out  into  four  parts:  contribution  to  the  equation  of  motion,  contribution  to  the 
boundary  condition  along  the  x  -  axis,  contribution  to  the  boundary  condition  along  the  y-axis, 
and  contribution  to  the  comer  boundary  condition.  These  expressions  are  represented  as: 

h 

fbUdt  =  ul  +  u2  +  u3  +  u4  (48) 

h 


Contribution  to  the  Equation  of  Motion  can  be  stated  as: 

h  bp.  a 

Ul  =  f  f  f  +  My’yydw0  +  2 M^dwjdxdydt 

q-fc/20 


(49) 


Contribution  to  the  Boundary  Condition  Along  the  X  -  Axis  can  be  stated  as: 


h  b/2 


U2  =  f  {  ~  Mx’x6w0  -  2Mxy’y6wO^ 


f,  -6/2 


(50) 


Contribution  to  the  Boundary  Condition  Along  the  Y  -  Axis  can  be  stated  as: 


U3  =  Jf WyM’y  ~  My’ydw0  ~  2Mxy’x5wO  I  ^ 

~bl  2 


(51) 


Contribution  to  the  Boundary  Condition  At  the  Comer  of  the  Plate  can  be  stated  as: 


-  a  6/2 

U4  =  f 2M  5w0 1  !  dt 

J  0-6/2 


(52) 


For  this  analysis,  the  structural  eigenvalue  problem  associated  with  harmonic  solutions 
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for  the  undamped  free  vibration  case  are  formulated.  Once  that  solution  is  found,  the  damped 
disturbed  vibration  case  will  be  formulated  based  upon  the  solution  to  the  above  problem.  The 
work  done  by  external  disturbances  can  be  defined  in  terms  of  conservative  torques  and  will 
show  up  in  the  boundary  conditions  for  the  damped  disturbance  vibration  case.  After  the 
equations  of  motion  are  solved  in  the  undamped  free  vibration  case  for  their  spatial 
relationships,  the  external  disturbances  and  the  control  force  are  adjoined  to  them  by 
LaGrange's  Principle.  Then  the  solutions  for  the  damped  disturbed  vibration  case  are  found. 

The  equations  of  motion,  initial  conditions,  and  boundary  conditions  can  now  be 
derived  by  substituting  the  appropriate  equations  from  the  variation  of  the  kinetic  and  potential 
energies  into  Hamilton's  Principle.  The  equations  are: 

Temporal  Boundary  Conditions  are: 


/  JM  .  l 

hKx^Wa\\dy  +  fe)> y6wol  I  dx 

-bp  *1°  0  h-b/2 


f  f  -  f3%’xx  ~  hKyy^Wo\  1^  =  0 


The  X-Axis  Boundary  Conditions  are: 


f  f  [-  Mx8w0,x  +  (Mx,g  +  2 M^y  -  hwQ,x)bw^ \dydt  =  0 
h-m  0 
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The  Y-Axis  Boundary  Conditions  are: 


//["  MybWpy  +  (My,y  +  |  dxA  =  0 


The  Comer  Boundary  Conditions  are: 


-  am  2 

-  f2M  6Wq|  |  =  0 

*  J  n-hn 


If  the  initial  and  final  conditions  of  the  problem  are  completely  specified,  then  the  variations 
of  the  degrees  of  freedom  over  the  interval  t,  to  t2  can  be  set  to  zero.  The  variations  of  the 
degree  of  freedom,  Sw0,  is  completely  arbitrary  and  is  not  necessarily  equal  to  zero.  Therefore 
the  quantities  multiplied  by  this  variation  must  be  set  equal  to  zero  in  order  to  satisfy 
Hamilton's  Principle.  The  partial  differential  equation  is  defined  as: 

Mx’xx  +  My’yy  +  2Mxy’xy  ~  73*W  "  hKyy  +  71^0  =  0  <58> 


No  assumptions  will  be  made  about  the  material  properties  until  the  equations  for  the 
piezoceramic  material  are  incorporated.  Therefore,  the  natural  and  geometric  boundary 
conditions  which  are  derived  from  Hamilton's  Principle  are  along  the  X  -  Axis: 


»*Wl  "  0 

0 


(Mx>x  +  2Mxy’y  ~  hKx^WoH  =  0  > 
J  0 
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along  the  Y  -  Axis: 


bp. 

(M)8w0  |  =  0 

J  -bp 


bp 


(My,y  +  2M^x  -  I3wQ,y) 8Wo]  I  =  0  , 


xy*x 


-bj  2 


(61) 

(62) 


and  at  the  comers: 

a  bp 

2 MJw0|  |  =  0  (63) 

^  0-b/2 

In  order  to  complete  the  development,  the  strain-displacement  expressions  and  the 
moment  resultants  need  to  be  substituted  into  the  above  equations.  Also,  the  material  constant 
I3  is  much  smaller  than  the  other  constants  in  the  equations  of  motion  and  boundary  conditions 
by  at  least  four  orders  of  magnitude.  Therefore,  this  constant  can  be  neglected  in  the  above 
equations.  Finally,  in  order  to  simulate  a  solar  array  structure  attached  to  a  satellite,  the 
boundary  condition  which  approximates  the  physical  conditions  is  a  fixed-free  condition  in  the 
x-axis  and  free-free  condition  in  the  y-axis.  This  is  a  satisfactory  boundary  condition  in  order 
to  control  the  relative  motion  of  the  solar  array  with  respect  to  the  satellite  bus.  Both  the 
displacement  w0  and  the  slope  of  the  displacement  w0,x  are  specified  at  the  fixed  end  at  (x=0), 
and- the  variation  of  these  variables  are  also  set  equal  to  zero.  For  the  free  ends  in  the  x-axis 
and  y-axis,  the  displacement  and  slope  cannot  be  specified  so  they  are  not  necessarily  zero. 
The  shear  and  moments  are  zero.  These  equations  can  now  be  rewritten  into  their  final  form 
in  terms  of  the  unknown  variables,  material  constants  and  proper  boundary  conditions.  The 
partial  differential  equation  is  defined  as: 
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^ll'^O’xxxt  +  2  (Pl2  +  2  +  4^l6W0’xxx, y 

+  D22WQ’yyyy  +  4D26W0^y  +  *1*0  =  0 

The  natural  and  geometric  boundary  conditions  become  along  the  X  -  Axis, 
at  X  =  0: 

6^(0)  =  0  "*  >W°)  =  0 

6w0(0)  =  0  -  wQ(0)  =  0 

at  X  =  a: 

Dnwo  .„(«)  +  DnWVyy(a)  +  2D16W0’Ja)  =  0 


16  O’Ay' 


DUW0’xJa)  +  (DU  +  4D66)WVxyy(a) 
+  4Dl6WQ^(a)  +  2D26W0’yyy(a)  =  0 


and  along  the  Y  -  Axis, 


at  Y  =  -b/2: 


Dl2W0’J-bl2)  +  D22W0’yy(-bl2)  +  2D26W^bl2)  =  0 

(£>12  +  4D26)w0,^(-i/2)  +  D22w0>m(-bl  2) 

+  2D16W0’xzx(~b/2)  +  4D66W0’*yy(~b/2)  =  0 


at  Y  =  b/2: 
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(71) 


A  ctive  Control  Derivation 

The  equations  of  motion  for  the  composite  plate  and  the  equations  governing  the 
active  control  elements  will  now  be  combined  to  form  the  complete  set  of  equations  that 
describe  the  structure.  Figure  1 1  outlines  how  the  active  control  elements  are  incorporated 
into  the  structural  equations  of  motion.  The  structural  equations  of  motion  are  combined  with 
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the  active  control  system  by  adjoining  the  equations  of  motions  with  the  sensor  and  actuator 
dynamic  model  equations.  The  performance  index  and  the  constraint  equations  are  formulated 
when  the  closed  loop  control  equations  are  formulated  from  the  combined  state  space  model. 

The  active  control  elements,  piezoceramic  actuators  and  sensors,  will  be  placed 
opposite  one  another  through  the  thickness  of  the  plate  and  will  be  wired  such  that  for  a  given 
voltage,  the  actuator  will  either  expand  or  contract  producing  a  bending  moment  at  the  actuator 
location.  They  will  primarily  be  used  to  control  the  bending  and  coupled  torsion  -  bending 
modes  of  vibration.  The  equations  for  the  bending  force  produced  from  the  piezoceramic  are 
now  stated.  (See  [8],  [39]  -  [45].)  These  equations  were  originally  derived  for  a  beam.  They 
are  applicable  in  this  case  due  to  the  dimensions  and  directions  of  the  capacitance  of  the 
actuators  and  sensors.  When  an  electrical  charge  is  applied  through  the  piezoceramics,  they 
will  either  expand  or  contract.  The  piezoelectric  coefficient  is  significantly  greater  along  the 
length  of  the  actuators  and  sensors  than  the  width.  Therefore,  they  will  produce  predominately 
a  pure  bending  moment  along  the  length  of  the  actuator.  These  equations  will  be  derived  for 
the  actuators  attached  on  the  top  and  sensors  attached  to  the  bottom  of  a  structure.  (See  figure 
12.)  Some  assumptions  made  in  deriving  these  expressions  are: 

1 .  A  perfect  bond  between  the  actuator  and  sensor  and  the  structure. 

2.  The  actuator's  and  sensor's  mass  do  not  contribute  to  the  structure. 

3.  The  capacitances  of  the  sensor's  leads  are  negligible. 

Considering  a  structure  with  no  external  moments  and  actuators  mounted  either  on  the  top  or 
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Figure  11  Combining  Equations  Of  Motion  with  Control  Theory 


ATTACHMENT  OF  SENSOR  AND 
ACTUATOR  TO  COMPOSITE 


z 


Figure  12  Actuator  and  Sensor  Attached  to  Structure 
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bottom,  the  expression  for  the  displacement  caused  by  these  actuators  can  be  expressed  as 
(see  [39]  -  [45]): 

w  =  — (77) 
El  2  - 


where  KD  is  the  piezoceramic  actuator  coefficient  and  VD  is  an  applied  DC  voltage.  This  is 
similar  to  the  expression  when  the  structure  has  no  actuators  present  and  is  acted  upon  by  a 
point  moment,  M0: 


w 


*v2 

2EI 


(78) 


Therefore,  the  piezoelectric  actuators  produce  an  equivalent  point  or  concentrated  moment 
proportional  to  the  applied  voltage.  As  stated  previously,  the  above  derivation  was  conducted 
for  beams.  Since  the  actuators  and  sensors  only  cover  a  small  portion  of  the  length  and  width 
of  the  plate,  the  moment  produced  and  the  strain  sensed  will  be  at  the  attachment  point.  This 
attachment  point  is  assumed  to  be  located  at  the  center  of  both  the  piezoelectric  actuator  and 
sensor. 

A  similar  expression  is  stated  for  the  piezoelectric  sensors  which  work  opposite  to  that 
of  the  actuators.  Instead  of  producing  a  strain  from  an  applied  voltage,  the  sensors  produce  a 
current  from  an  applied  strain  rate.  The  results  of  these  sensor  equations  found  in  references 
[39]  -  [45]  are: 


*  r?  ' 

VAt)  =  Ks  f  z—(x,t)dx  ,  k  =  1,2,  ...,  N 
*  /  dxz 


(79) 
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After  integrating,  this  equation  becomes: 


V  (t)  =  | 

k  dX  X; 


k  =  1,  2,  IV 


(50) 


where: 

K,  =  (S ■ 

and 

v0,k 

= 

Measured  voltage  from  the  kth  sensor 

wd 

= 

sensor  width 

h/2 

= 

one  half  plate  thickness 

N 

= 

number  of  sensor  and  actuator  pairs 

d31 

= 

piezoelectric  constant 

cEu 

= 

elastic  modulus  of  the  sensor  material 

*k-l 

= 

beginning  point  of  kth  sensor 

*k 

= 

endpoint  of  the  kth  sensor 

Rf 

= 

resistor  associated  with  an  ideal  operational  amplifier  used  to  acquire 

the  signal 

In  this  research,  N=2,  two  actuator  and  sensor  pairs  will  be  used  and  their  locations  and 
control  gains  will  be  optimization  parameters. 

Discretization  of  the  Equations  of  Motion 

The  equations  of  motion  for  the  entire  plate  including  the  active  control  system  will 
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now  be  reduced  to  ordinary  differential  equations  with  respect  to  time  to  formulate  the  closed 
loop  state  space  equations  based  upon  the  control  theory  used.  From  these  closed  loop  state 
space  equations,  the  performance  index  and  constraint  equations  are  developed.  This  can  be 
done  by  determining  the  approximate  solution  to  the  spatial  dependence  and  integrating  this 
solution  over  the  volume  of  the  plate  leaving  only  an  ordinary  differential  equation  with 
respect  to  time.  An  assumption  concerning  the  above  equations  must  be  made  to  simplify  the 
expressions.  This  assumption  is  that  the  spatial  dependence  can  be  separated  in  the  x  and  y 
directions  and  from  the  time  dependence  terms  in  the  equations  of  motion.  This  assumption  is 
hard  to  prove  since  it  assumes  that  synchronous  motion  occurs  from  which  orthogonality  is 
shown.  Taking  into  account  these  assumptions,  the  x  and  y  spatial  coordinates  and  the  time 
dependence  are  separated  and  an  approximate  method  such  as  the  assumed-modes,  Rayleigh- 
Ritz  method  (see  [33],  [35],  [36],  [46]  -  [49])  or  Galerkin  (see  [33],  [35],  [36],  [46],  [47],  [49]) 
is  used  to  solve  the  equations  of  motion.  A  typical  Rayleigh-Ritz  approimate  solution  for  this 
application  is: 

N  M 

wQ(x,y,t)  =  £  £  (82) 

i  =  1 }  =  1 

Notice  the  spatial  and  temporal  dependencies  are  separated.  This  approximate  spatial 
relationship  or  admissible  function  is  assumed  to  satisfy  only  the  geometric  boundary 
conditions.  For  this  research  with  the  above  boundary  conditions  (eqs  65-76),  fixed-free  along 
the  X-axis  and  free-free  along  the  Y-axis,  geometric  boundary  conditions  are  present  only 
where  the  plate  is  fixed  in  the  X-axis  direction.  The  spatial  properties  or  geometric  boundary 
conditions  for  the  X-axis  are  satisfied  by  the  admissible  function  [33],  [49],  [50]: 


41 


XJM  =  1  -  cos[^H] 
2  a 


m  =  1,2, 3, 4,. 


(83) 


For  the  Y-axis,  the  function  used  is  a  comparison  function  for  the  flexible  motion  of  a  beam 
made  out  of  an  isotropic  material  with  free-free  boundary  conditions  [33],  [50],  This  can  be 
derived  from  assuming  motion  is  the  sum  of  four  trigonometric  functions:  cos,  sin,  cosh,  and 
sinh. 

Y  =  c.sin—  +  c,  cos—  +  c,siith—  +  c.cosh—  (84) 

1  b  2  b  3  b  4  b 


where 

Cp  c2,  c3,  c4  constants  derived  from  the  boundary  conditions 

X  solution  to  the  characteristic  equation 


Substituting  the  above  expression  into  the  boundary  conditions  for  free-free  which  are: 


— (0) 

ay2 

0(i) 

dy 


&Y 
dy 3 

&Y 

dy3 


(0)  =  0 
(b)  =  0 


(85) 


yields  the  following  four  expressions: 
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(86) 


X2 

—S~C1  +  C4>  =  0  C2  =  C4 

bz 


-U-c,  ♦  c3)  =  0  - 

tr 

C1  =  C3 

(87) 

(-CjSinA  -  c2 cosA  +  cqsmliA 

+  c2coshA)  =  0 

(88) 

(-CjCosA  +  c2sinA  +  CjCoshA 

+  c2siuhA)  =  0 

(89) 

Solving  equations  (88,  89),  assuming  X  doesn't  equal  0  which  would  be  the  trivial  solution, 
one  obtains: 

cosAcoshA  =  1  X  =  4.73,  7.8532,  10.996,  14.137,....  (90) 

and 

[coshA  -  cosA] (91) 
1  [sinhA  -  sinA]  2 


Using  the  above  equations,  the  expression  for  the  flexible  motion  of  a  beam  made  out  of  an 
isotropic  material  with  free-free  boundary  conditions  is: 

Y  =  cosh[-^]  +  cos[^] 


(92) 


(cosh[A]  -  cos[A]) 
(sinh[A]  -  sin[A]) 


(sinh[AZ]  + 


sin[^]) 

b 


However,  additional  terms  are  needed  since  the  assumed  comparison  function  modeled 
only  the  flexible  motion  of  the  beam  and  the  rigid  body  motion  of  the  beam  in  the  Y-axis  also 
needs  to  be  taken  into  account.  These  additional  terms  capture  the  rigid  body  motion  of  the 
plate  in  the  Y-axis  while  the  plate  in  the  X-axis  is  undergoing  flexible  body  motion.  Also  in 
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order  to  simplify  the  equations,  the  coordinate  system  for  the  Y-axis  was  shifted  from  the 
center  of  the  plate  to  one  end.  This  shift  occurred  over  a  distance  of  b/2  or  half  the  width  of 
the  plate.  (See  figure  13  for  the  new  coordinate  system.) 


r,  *  i 

r2. 1  -  2 

2  b 

Yn(y)  =  cosh[— ]  +  cos[^rf] 
b  b 

_  (cosh[A(f,_2)]  -  cos[A(,_2)])  A(n,2)y  + 

(sinh[A(i|_2)]  -  sin[A(l|_2)])  b  b 


]) 


n  =  3,4, 


(93) 


where  An  is  a  solution  to  the  characteristic  equation  given  in  equation  (90).  The  approximate 
periodic  solution  is  a  combination  of  both  solutions  and  \|/(t)  =  sin[pt]: 


m  n 

W0 (x,y,t)  =  £  £  A^XYsmlpt] 
i  =  1  j  =  1 


(94) 


using  a  finite  number  of  terms  for  X  and  Y  (in  this  case  2  terms  each): 


w0(x,y,t)  =  [An(l  -  cos^)  +  X12(l  -  cos^)(l  - 

2  a  2  a  b 

+  A2y(  1  -  cos—)  +  X22(l  -  cos— )(1  -  —  )]smpt 
a  a  b 


(95) 


and: 

p  =  Modal  frequency 

Following  the  Rayleigh-Ritz  approach,  this  approximation  is  then  substituted  into  the  following 
expression  for  kinetic  and  potential  energy  from  which  the  equations  of  motion  for  the 
undamped,  unforced  plate  were  derived  to  determine  the  frequencies  and  mode  shapes  of  the 
plate: 
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■5 


(96) 


-  r_ 


=  o 


where: 


max 

t 


b  a 

<i/A 
0  0 


73w0^ 


/3w0)J 


^wjjdrafy) 


(97) 


h 

b  a  2 

^max  =  “T  If  f  (°A  +  °,S  +  ^xyVxy^^y 

0  0  Jt 
2 


(98) 


b  a 

max  r  r  1 


/  /  -HM*W0,xx  +  +  2MxyWoJ^ 


f  J  j  2 
0  0 


These  are  used  to  solve  for  the  constants,  Ay,  and  the  frequencies  of  the  approximate  solution. 
After  the  integrations  with  respect  to  x,  y,  and  z  are  completed  and  the  numerical  values  are 
substituted  in  for  the  constants  found  in  figure  10,  the  expressions  for  the  potential  and  kinetic 
energies  become: 

=p2?[MA2n  +  .184  +  2.8AuA21  +  3.64  +  l-^J] 

Vma  =  -00044 A2nDn  *  .00015A*DU  *  .0015 AnAnDn 

+  ,007L4221Du  +  ,0005A12A22Du  +  ,0023A222DU 
,0QS6AnAl2D16  +  .023 Al2A21D16  -  .023A11A22D16 

+  A0A  ftDfa  +  ■^■^1^-22^66  +  -^^22^66 
The  derivative  of  the  above  equation  is  taken  with  respect  to  the  constants  Ay  and  each  is  set 


(99) 


(100) 
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equal  to  zero.  This  is  used  [36]  to  calculate  the  constants  Ay  from  which  the  eigenvector  can 
be  determined.  This  forms  four  separate  equations  in  terms  of  the  unknown  Ay  and 
frequencies.  These  four  equations  are  combined  into  a  matrix/vector  equation  in  order  to  solve 
for  the  constants  and  frequencies.  The  constants  and  frequencies  are  then  found  by  solving  the 
equations  simultaneously. 

The  method  that  is  used  in  this  research  is  evaluated  by  differentiating  each  term 
separately  and  formulating  the  following  equation: 


dTmax 

8An 

MU 

3^ 

34ax 

U12 

3^ 

►  —  - 

3^12 

3T-J 

3A2i 

3^21 

3^max 

8T 

max 

8A22 

3A22 

(101) 


Each  term  can  be  represented  as  a  matrix  multiplication  using  the  symbols  K  and  M  for  the 
stiffness  and  mass  matrices  respectively. 


dV 


dA 

dV 


11 


dA 

0K 


12 


dA. 

3F 


21 


dA 


22 


=  m 


21 


22 


(102) 
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BT 


BA 

BT 


n 


aA 

BT 


12 


at 

BT 


21 


a4 


22 


pi 


=  p2[M] 


11 


22 


(703) 


where 


.00088DU  -.0086Di6 

.0015DU 

-,023D16 

-.0086D16  .00029Du+.21D66 

,023D16 

.0005£>U+.35D£>& 

,0015DU  ,023D16 

.014Dn 

0 

-,023£»16  .000 SDn+,351)^ 

0 

,0047£>u+.82D66 

(104) 


and 


[M]  =  p 


1.1  0 
0  .36 

2.76  0 
0  .92 


2.76 

0 

0 

.92 

7.2 

0 

0 

2.4 

(205) 


Substituting  these  expressions  into  the  above  equation,  one  obtains  the  following  equation: 


M 

[K\ 

^12 

^21 

•  =  p2[M]‘ 

An 

^21 

^22 

422 

(106) 


The  frequencies  and  constants  A-  are  found  by  solving  the  generalized  eigenvalue  problem  for 
the  above  equation.  The  frequencies  squared  are  the  eigenvalues  and  the  constants  Ay  are  the 
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eigenvectors. 

As  a  check  to  see  if  the  selected  approximate  solution  models  the  plate  sufficiently,  the 
numerical  frequencies  and  mode  shapes  were  found  by  using  typical  realistic  material  constants 
and  various  ply  layups.  For  this  check,  typical  composite  parameters  were  substituted  for  the 
composite  properties.  (See  Table  1.) 


TABLE  1  TYPICAL  COMPOSITE  PROPERTIES 


Composite  Parameter 

Value 

E, 

137.25  GPa 

e2 

11.42  Gpa 

G  jo 

6.39  GPa 

V,2 

0.3 

V21 

0.025 

p 

.043 

These  results  for  the  four-mode  approximation  (2  "x"  terms  and  2  "y"  terms)  also  compare,  at 
least  in  first  bending,  to  the  one-mode  approximation.  The  frequencies  for  the  one-mode 
approximation  are  now  calculated  for  the  following  four  test  case  ply  layups:  1.  [0/90/±45/0]s, 
2.  [0/±452]s,  3.  [±455],  4.  [+305]s  and  shown  in  table  2.  They  differ  by  between  1-3%  for  all 
four  cases.  Since  the  Rayleigh-Ritz  solution  is  an  upper  bound,  the  frequencies  should 
decrease  for  the  more  accurate  four-mode  approximation. 
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TABLE  2  FREQUENCIES  FOR  THE  TEST  CASES  FOR  THE  ONE  MODE 

APPROXIMATION 


Case  # 

Frequency  (Hz) 

1 

1.3857 

2 

1.03544 

3 

1.46278 

4 

1.0116 

These  are  the  frequencies  for  the  first  bending  mode  of  the  structure.  The  frequencies  and 
associated  mode  shapes  are  now  calculated  for  the  four-mode  approximate  solution  for  the 
above  test  cases  and  shown  in  table  3. 


TABLE  3  FREQUENCIES  AND  MODE  SHAPES  FOR  THE  TEST  CASES  FOR  THE 

FOUR  MODE  APPROXIMATION 


Case  # 

Frequency  (Hz) 

Mode  Shape 

1 

1.3416 

1st  bending 

11.207 

1st  torsion 

12.45 

2nd  bending 

49.11 

2nd  torsion 

2 

1.00831 

1st  bending 

9.2879 

1st  torsion 

18.745 

2nd  bending 

81.138 

2nd  torsion 

3 

1.41583 

1st  bending 

12.9074 

1st  torsion 
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14.7732 

2nd  bending 

63.697 

2nd  torsion 

1.00395 

1st  bending 

7.535 

1st  torsion 

17.67 

2nd  bending 

73.43 

2nd  torsion 

The  mode  shapes  are  what  one  would  have  expected  for  a  plate  geometry.  However,  it  is 
difficult  to  prove  this,  because  there  is  no  known  solution  for  this  problem  and  no  one  has 
modelled  this  particular  plate  before.  Other  people  have  modelled  plates  with  fixed-free  end 
conditions  before,  but  one  was  not  found  for  these  particular  ply  layups  and  plate  dimensions 
[9],  [33],  [50],  [51].  For  a  plate  structure  which  is  very  long  in  one  direction  and  very  thin, 
one  would  expect  the  low  frequency  mode  shapes  to  be  in  the  dominant  direction.  Also  with 
the  purposeful  ply  angle  selection  to  enhance  the  bending/torsion  coupling,  one  would  also 
expect  the  combined  bending/torsional  mode  shape  to  appear  in  the  first  few  modes  of 
vibration.  Furthermore,  the  frequency  changes  among  the  different  test  cases  (ply  angles) 
make  sense.  The  more  rigid  the  plate  is  in  torsion,  the  more  the  frequency  for  that  particular 
mode  will  increase.  Also  if  there  are  0  degree  plies  in  the  layup,  the  frequencies  for  the 
bending  modes  will  increase  due  to  the  increased  longitudinal  stiffness  of  the  plate.  Based 
upon  the  calculated  frequencies  and  mode  shapes,  it  is  assumed  that  the  approximate  solution 
models  the  plate  sufficiently  for  this  research.  As  a  final  check,  the  frequencies  of  a  nine¬ 
mode  approximation  (3  "x"  terms  and  3  "y"  terms)  are  solved  numerically  for  case  number  1. 
These  results  are  compared  with  the  one  and  four  mode  approximations  for  case  1.  The 
frequencies  for  all  modes  are  presented  in  table  4. 


51 


TABLE  4  TEST  CASE  1  FREQUENCIES  FOR  THE  NINE  MODE  APPROXIMATION 


Mode  Number 

Frequency  (Hz) 

1 

1.31765 

2 

8.736 

3 

10.812 

4 

34.25 

5 

49.838 

6 

102.59 

7 

320.36 

8 

327.5 

9 

365.94 

Table  5  shows  the  comparison  of  the  first  four  modal  frequencies  for  the  one,  four,  and  nine 
mode  approximations. 
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TABLE  5  TEST  CASE  1  FREQUENCIES  FOR  THE  ONE,  FOUR,  AND  NINE  MODE 

APPROXIMATION 


Mode  Shape 

One  Mode 

Four  Mode 

Nine  Mode 

Percent  Difference  4  vs  9  Mode 

1st  Bending 

1.3857 

1.3416 

1.31765 

1.8% 

2nd  Bending 

X 

11.207 

8.736 

22% 

1st  Torsion 

X 

12.45 

10.812 

13.2% 

X 

49.11 

34.25 

30% 

As  one  can  see,  the  frequencies  for  the  nine  mode  approximation  are  lower  than  both  the  one 
and  four  mode  approximation.  The  frequencies  for  the  four  mode  approximation  vary  from 
the  nine  mode  approximation  by  2-30%.  Clearly,  the  nine  mode  approximation  would  be  a 
better  solution  to  use,  but  a  symbolic  solution  was  never  calculated  due  to  the  limitation  of  the 
computer  workstations.  It  is  assumed  that  the  four  mode  approximation  is  sufficient  for  this 
research.  Also  in  Chapter  VI,  this  approximation  will  be  compared  to  experimental  results  and 
the  referenced  paper  by  Hwang  [25]  which  used  a  finite  element  method.  These  comparisons 
will  show  that  this  approximate  solution  models  the  plate  sufficiently. 

The  Rayleigh-Ritz  solution  using  generalized  time  functions,  versus  periodic 
expressions  for  time,  for  the  modal  amplitudes  can  now  be  formulated  using  the  approximate 
solution.  The  dimensions  of  the  plate  are  substituted  into  the  appropriate  equations  prior  to  the 
formulation.  This  solution  is: 

w0(W)  =  wfcjOjil^O)  +  w(x,y)2^2(t)  +  W(x,y)3^3(t)  +  w(x,y)4i]t4(0  (107) 
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where: 


MxA  -  Au&r,  *  AvJCJi  *  *  A22tX2Y2  (108) 

and: 

X  =  (1  -  cos^)  m  =  1,2  (109) 

m  48 

Yx  =  1  (UO) 

Y2  =  (1  -  2)  (Hi) 

2 

In  the  above  equations,  the  numerical  plate  dimensions  were  used  versus  the  symbols,  a  and  b. 
The  constants  Ay  vary  for  each  spatial  solution,  w(x,y)j.  These  assumed  modes,  w0(x,y,t),  can 
not  be  substituted  directly  back  into  the  equations  of  motion  because  not  all  of  the  boundary 
conditions  are  satisfied.  Only  the  geometric  boundary  conditions  are  satisfied.  (Note  that  an 
admissible  function  versus  a  comparison  function  was  used  in  the  Rayleigh-Ritz  approximation 
so  only  the  geometric  boundary  conditions  should  be  satisfied.)  If  we  use  only  the  equations 
of  motion,  then  we  completely  ignore  the  natural  boundary  conditions.  The  assumed  modes 
are  then  substituted  into  the  expressions  for  potential  and  kinetic  energies  and  the  spatial 
dependence  is  integrated  out.  To  generate  the  new  equations  of  motion  we  can  use  either 
Hamiltion's  Principle  with  the  \|/;(t)'s  as  independent  variables  or  the  Lagrange  equation,  since 
we  are  dealing  with  conservative  forces  and  motions  in  nature.  For  ease  of  use  to  generate  the 
new  equations  of  motion  and  the  ability  to  incorporate  external  forces,  the  Lagrange  equation 
is  chosen.  The  Lagrange  equation  will  generate  four  equations  of  motion  for  each  of  the 
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independent  variables. 


d(  dT  " 


dT  +  dV 

#,(0 


i  =  1 ,2,3,4 


(112) 


The  complete  equation  of  motion  for  the  plate  is  formed  by  combining  the  four  equations  of 
motion  into  a  matrix  equation: 

[M]$(t)  +  [jqi|r(f)  =  Q  (]J3) 


where  M  and  K  are  given  in  equation  (102)  and  (103) 

i|t(0  =  {i|r(f),-}  .  (II4) 

and  vjr(t)  represents  the  time  dependent  modal  amplitudes  of  the  approximated  solution.  The 
next  step  is  transforming  this  set  of  ordinary  differential  equations  into  a  state  space 
representation.  This  is  done  in  the  next  section  in  conjunction  with  adjoining  the  control 
system  dynamics  to  the  flexible  body  equations  of  motion. 

Control  Theoiy  and  Optimization 

In  this  section  the  equations  of  motion  from  the  previous  section  are  transformed  from 
a  second-order  ordinary  differential  equation  into  a  state  space  format  that  is  a  first  order 
ordinary  differential  equation  and  then  incorporated  into  a  control  theory  for  the  active  control 
system.  A  performance  index  and  constraint  equations  which  formulate  the  optimization 
problem  are  then  determined  by  using  the  state  space  equations. 

Control  Theory.  The  control  theory  used  in  this  research  is  output  feedback  from 
collocated  sensors  and  actuators.  Collocated  is  defined  as  occupying  the  same  location  on  the 
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plate  along  the  X-axis  and  Y-axis.  The  sensors  and  actuators  can  be  located  on  opposite  sides 
of  the  plate  and  still  be  collocated.  Output  feedback  is  defined  as  taking  the  signal  obtained 
by  a  sensor,  multiplying  the  signal  by  a  gain  and  transmitting  it  to  a  control  actuator  to 
produce  a  force  which  could  be  a  moment.  This  sensor  senses  the  strain  rate  of  the  plate.  An 
active  control  electronic  filter  is  sometimes  needed  to  process  the  output  signal  that  comes 
from  the  sensors  and  is  converted  into  a  high  voltage  control  signal  that  is  provided  to  the 
actuators.  This  filter  is  used  only  to  eliminate  any  high  frequency  noise  which  is  also  detected 
by  the  sensor.  A  standard  active  control  feedback  filter  used  with  piezoceramic  sensors  and 
actuators  consists  of  a  preamplifier,  a  bandpass  filter,  a  phase  shifting  circuit,  and  a  high 
voltage  output  amplifier.  A  filter  was  needed  and  the  center  frequency  was  set  to  the  torsional 
frequency  of  the  plate.  From  the  definition  of  direct  output  feedback,  one  sensor  controls  only 
one  actuator. 

The  actuator  moments  and  disturbance  torques  are  added  to  the  equations  of  motion 
now  by  equating  the  individual  equations  of  motion  resulting  from  LaGrange's  equation  to  a 
work  term. 


d(  dT  )  _  dT  +  dV 
dt{  3t[r  (.(r)  J  3t|r,.(f)  3^.(0 


(U5) 


where: 

N 

<?,  -  £  ft 

k  =  1 

The  force,  Fk,  in  the  above  equation  is  the  point  actuator  moment  applied  anywhere  on  the 

— ?■ 

structure  where  the  actuator  sensor  pair  is  located.  The  force,  Fj,  in  the  above  equation  is  the 
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disturbance  torque  applied  to  the  structure.  Since  the  force  is  treated  as  a  moment,  the  vector 
"q  is  the  vector  whose  magnitude  is  the  angle  over  which  the  moment  acts.  For  the  actuators, 
the  vector!^  is  the  vector  whose  magnitude  is  the  angle  or  X-axis  rotation  of  the  plate  where 
the  actuator  is  located  at  x  =  xAk  and  y  =  yA  k.  This  is  modelled  as: 


i  act 


=  m. 


dw 

dx 


dijt,. 


•*  =  y*k 


(117) 


K„Vn 


(118) 


The  disturbance  torque  is  applied  to  the  free  end  of  the  plate  producing  a  twisting  of  the  plate. 
Again,  since  the  force,  Fj(  is  a  torque,  the  vector  Tj  is  the  vector  whose  magnitude  is  the  angle 
over  which  the  moment  acts.  This  is  modelled  as: 


disturbance 


=  *D  » 


y  =  xd 


(119) 


where  d  is  the  disturbance  torque  amplitude  and  xD  =  a  and  yD  =  b/2.  The  vector  p,  in  this 
case,  is  the  vector  whose  magnitude  is  the  angle  or  Y-axis  rotation  of  the  plate. 

Before  reducing  the  ordinary  differential  equations  to  state  space  equations  the  above 
actuator  and  disturbance  models  are  incorporated  into  the  model. 


[W  +  im  = 


{  dw^ 

N  u 

Y''  K  V  _ 

u*, 

i 

+  /UJ 

1 

k  =  t 

5v(f 

U  =  Xjtk 

y  =  3t|r 

=xD>y  =y0 

(120) 


The  sensor  equation  makes  up  the  last  part  of  the  state  space  transformation.  This  equation  is: 
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V-S »  =  vf  H 


where  xk  is  the  end  point  and  xk_,  is  the  beginning  point  of  the  sensor  corresponding  to  the 
actuator  location.  The  second  order  equation  of  motion  is  rewritten  as: 

$  =  -[MT1^  +  [Ml~\b]  [FbJ 


_i[ — )  3  d 

L  a,lr  hr  =  Xp  ,  y  =  yDi 


where 


M  =  W 


is  a  vector  of  the  control  voltages  and  the  iklh  element  of  [b]  is: 


The  state  variable  q  is  related  to  the  physical  variables  with: 
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The  second  order  differential  equation  can  now  be  rewritten  in  first  order  state  space  notation 
in  which  the  vector  q  represents  the  assumed  modal  displacements  and  first  derivatives.  The 
first  order  state  space  notation  is  given  by: 


q  =  [A]q  +  [B]w  +  [H]d 


(127) 


y  =  [C\q 


(128) 


where: 


-[MP[£]  0. 


(129) 


B  = 


0 


(130) 


0 


H  = 


=  xd  >  y  =  yD 


(131) 


C  =  [0  [c]} 


(132) 


where  the  kith  element  of  [c]  is: 


■  V 


[ax. 

d$i 


xk-i 


(133) 
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and: 


[A]  =  Plant  matrix  determined  from  the  stiffness  matrix  and  the  inverse  of 

the  mass  matrix 

[B]  =  Control  matrix  determined  from  the  actuator  coupling  coefficients, 

location  matrix  and  the  inverse  of  the  mass  matrix 
u  =  [VDk],  Control  input  vector  which  is  the  actuator  voltages  at  specific 

locations 

y  =  Sensor  voltage  signal  (strain  rate  of  the  plate)  collocated  at  actuator 

position 

[C]  =  Sensor  matrix  determined  from  the  sensor  and  dynamic  coupling 

coefficients 

[H]  =  Disturbance  matrix  determined  from  the  disturbance  location  matrix 

and  the  inverse  of  the  mass  matrix 
d  =  Torsional  disturbance  amplitude  at  free  end  of  plate 

In  output  feedback,  the  control  vector  u  (voltages  sent  to  the  actuator)  is  proportional  to  the 
output  sensor  vector  y  (strain  rate  of  the  plate)  by  a  gain  multiplier. 

u  =  -[G]y  (134) 


where: 

[G]  =  constant  control  gains 

Substituting  equation  (128)  into  the  expression  for  y,  one  obtains: 
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k  =  ~[G\\c\q 


(135) 


Substituting  the  above  equation  into  the  state-space  relationship,  the  closed  loop  state-space 
equation  can  now  be  derived. 


q  =  [A]q  -  [B][G][C]q  +  [H]d  (136) 

q  =  [AJq  +  [H]d  (137) 

where: 

[AJ  =  [A]-[B](G\[Q  (138) 

Performance  Index  and  Sensitivity  Equations.  The  performance  index  or  objective 
function  is  a  scalar  measure  that  is  to  be  minimized  during  the  structural/control  optimization. 
The  primary  objective  of  the  active  damping  system  is  to  minimize  the  bending  and  torsional 
motion  of  the  plate  while  staying  within  the  actuator  voltage  constraints.  Previous  papers  that 
use  output  feedback  as  their  control  system  tend  to  minimize  the  norm  of  the  output  feedback 
gain  matrix  [54]  which  was  a  function  of  the  structural  and  control  parameters  or  the 
expectation  of  the  standard  Linear  Quadratic  Gaussian  plus  the  penalized  gain  matrix  [39]. 
This  research  will  be  similar  to  the  work  accomplished  by  Slater  [56]  in  which  he  states  "the 
structure-control  design  is  to  find  the  structural  parameters  and  the  control  law  to  minimize  a 
performance  index  while  satisfying  control  energy  and  displacement  constraints."  The 
structural  equations  of  motion  are  assumed  disturbed  by  a  torsional  zero  mean  white  Gaussian 
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noise  disturbance  with  an  intensity  of  D.  A  white  Gaussian  noise  was  chosen  as  the 
disturbance  to  optimize  the  performance  of  the  plate  to  a  wide  range  of  frequencies  versus  one 
set  of  frequencies  as  in  harmonic  motion. 

The  performance  index  or  objective  function  J  used  in  this  research  will  be  based  upon 
the  mean  square  motion  of  the  free  end  of  the  plate  wtip  (ps),  where  the  pj  are  the  structural  and 
control  parameters  to  be  optimized,  e.g.,  the  ply  angle,  location  of  the  actuators,  and 
magnitude  of  the  actuator  gains.  The  stiffness  of  the  plate,  which  reacts  against  the 
disturbance  force,  is  determined  from  the  lamina  stiffness  which  depends  on  the  ply  angle. 

Also  the  coupling  parameter  is  determined  directly  from  the  ply  angle.  Just  increasing  the  ply 
angle  does  not  cause  the  coupling  parameter  to  increase  as  shown  by  the  equations  for  the 
lamina  stiffness.  The  coupling  parameter  will  reach  a  relative  maximum  for  some  optimum  ply 
angle. 

The  determination  of  this  optimum  angle  is  shown  next.  This  will  show  that  there  is 
only  one  ply  angle  which  maximizes  the  torsionai/bending  coupling  coefficient,  Dl6.  Because 
the  Y-axis  approximation  used  in  this  research  is  linear,  the  D26  coefficient  is  not  taken  into 
account  when  calculating  the  equations  of  motion.  The  D26  coefficient  is  multiplied  by  the 
second  derivative  of  the  displacement  function,  w0(x,y,t),  with  respect  to  y  when  forming  the 
potential  energy  equation.  This  expression  is  zero  due  to  the  linear  Y-axis  displacement 
function.  (See  equations  110,111.)  The  equation  for  the  D16  coefficient  can  be  expressed  as: 

Die  =  EiQJl]  z2dz  (m 


where: 
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(140) 


Ql2@26 


-22 


and  [Qjj]  given  by  eqs  (14-19).  For  this  problem,  the  experimentally  determined  values,  with 
the  exception  of  v21  which  is  calculated,  used  for  the  material  properties  are: 

E,  =  156.45  GPa 
E2  =  9.95  GPa 
G12  =  5.99  GPa 
v12  =  .3056 
v21  =  .01925 

The  angle  which  optimizes  the  D16  coefficient  can  be  determined  by  plotting  this  expression 
for  magnitude  D,6  versus  0  or  by  taking  the  first  derivative  of  the  expression  for  D16  and 
setting  it  equal  to  zero  and  determining  0.  (See  Figure  14.)  The  angle  (expressed  in  radians) 
which  maximizes  the  D16  coefficient  is: 

0opt  =  .312 
or  expressed  in  degrees: 

0opt  =  18.0° 

The  objective  of  the  optimization  problem  in  this  research  was  to  find  the  feedback 
law -or  gain  G,  actuator  locations  and  ply  angles  which  minimize  J  subject  to  certain 
constraints  which  are  related  to  the  mean  square  control  energy.  The  mean  square  control 
energy  constraint  provides  an  upper  limit  on  the  amount  of  control  energy  expended  since  the 
piezoceramic  actuators  do  have  an  upper  limit.  This  optimization  problem  can  be  stated  as: 

Minimize: 

J  =  CE(wrwrT)  *  (1  -  QE(wbwTb)  (141) 
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D16  Coefficient 


MAXIMUM  D16  COEFFICIENT 


Figure  14  D16  Coefficient  versus  Ply  Angle 
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subject  to: 


E(uut)  <  [S] 


(142) 


where: 

G 

Parameter  to  vary  the  weight  on  a  torsional  versus  bending 

motion 

wr 

Difference  of  the  corner  displacement  vector  and  the  center 

displacement  vector  of  the  plate  tip 

wb 

Center  displacement  vector  of  the  plate  tip 

E 

Fixed  mean  square  control  energy 

E()  = 

Expectation  operator 

u  = 

Control  vector  {actuator  voltages} 

The  performance  index  and  inequality  constraints  can  be  converted  from  statistical 
notation  in  terms  of  the  expectation  operator  to  matrix  notation  using  the  Lyapunov  equation 
solution  of  the  covariance  matrix.  The  motion  of  the  tip  wr  and  wb  are  determined  from 
multiplying  the  location  matrix  by  the  state  response: 


=  [Cr]q 


wb  = 


(143) 


For  large  time  behavior  the  motion  can  be  represented  by  the  steady  state  covariance  of  the 
states: 
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E(qqT)  =  Q 


(144) 


The  steady  state  covariances  can  be  calculated  by  solving  a  Lyapunov  equation  [56],  Both  the 
performance  index  and  the  constraint  equations  can  be  rewritten  in  terms  of  matrix  equations 
due  to  the  relationship  between  the  mean  square  response  of  the  states  and  the  covariance 
matrix  Q.  The  mean  square  responses  of  the  states  will  be  the  diagonals  of  the  covariance 
matrix  Q.  Before  the  Lyapunov  equation  is  solved,  the  closed  loop  plant  matrix  is  calculated. 
The  closed  loop  plant  is  solved  by  substituting  the  output  and  control  expression  into  the 
equation  of  motion.  This  is  represented  as: 

[AJ  =  [A]  -  [B][G][q  (145) 

The  Lyapunov  equation  [56]: 

[AJ[Q]  *  [Q][AJT  +  [H][D][HiT  =  0  (m 

is  solved  for  the  covariance  matrix  Q,  where  D  is  the  intensity  of  the  white  Gaussian  noise. 
For  steady  state  optimization,  the  performance  index  and  actuator  energy  constraints  can  be 
replaced  by  their  respective  covariance  relationships.  The  performance  index  is  rewritten  in 
matrix  form  by  substituting  the  covariance  matrix  in  for  the  expected  value  of  the  state: 

E(wrwrT)  =  E(Crqq  TCj) 

=  C,E(qqT)C?  (147) 

=  CfiCj 
=  E(Cbqq  TCTb) 

=  CbE(qqT)CbT  (148) 

=  CbQCbT 
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The  constraint  equations  are  rewritten  in  matrix  form  by  substituting  the  covariance  matrix  in 
for  the  expected  value  of  the  state: 

E(u  uT)  =  E(GCqqTCTGT) 

=  GCE(qqT)CTGT  (149> 

=  GCQCtGt 


The  constraint  equations  form  a  two-by-two  diagonal  matrix  whose  elements  are  required  to  be 
less  than  or  equal  to  the  total  control  energy  the  actuators  can  produce. 


GCQC'G 


Tr'T  _ 


Cj  0 


0  c, 


(150) 


The  total  control  energy,  s,  is  also  a  two-by-two  diagonal  matrix  whose  elements  are  equal 
and  equate  to  the  allowable  control  energy  each  actuator  can  produce. 


'«!  0 
0  S2 


(151) 


The  performance  index  and  constraint  equation  can  be  finally  rewritten  as: 
minimize 

J  =  acrMC/  +  (1  -  0[CJ<2[C/  (152) 


subject  to 

cl<Zl 

c2<?2 


(153) 
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where: 


c 

cr 


Cb 

Q 


^2 

c 

G 


Parameter  to  vary  the  weight  on  a  torsional  versus  bending 
motion 

Location  matrix  for  the  difference  between  the  left  comer  and 
the  center  of  the  plate  tip 
Location  matrix  for  center  of  the  plate  tip 
Covariance  matrix  from  Lyapunov  equation  (152) 

Fixed  voltage  control  energy 

Fixed  mean  square  control  energy  for  actuator  one 

Fixed  mean  square  control  energy  for  actuator  two 

Location  matrix  for  sensor  output 

Actuator  gain  matrix 


The  optimality  conditions  or  sensitivity  equations  which  arise  from  the  optimization 
are  calculated  by  taking  the  partial  derivative  of  the  performance  index  and  the  constraint 
equations  with  respect  to  the  parameters  to  be  optimized.  These  parameters  for  this 
performance  index  are  the  structural  and  control  parameters  pj  which  include  the  actuator 
locations,  and  the  control  gains  G.  These  gradients  are  calculated  numerically  by  the  finite 
difference  method  by  an  optimization  program  called  DOT  [55],  The  search  gradients 
determine  which  direction  the  optimization  program  proceeds.  This  will  be  discussed  in  the 
next  chapter  on  optimization. 
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IV.  Optimization  Results 


A  FORTRAN  computer  optimization  program  called  DOT  [55]  was  used  to  solve  the 
optimization  problem  in  this  research  using  numerical  optimization.  Numerical  optimization 
determines  the  "best"  solution  using  a  logical,  systematic  decision-making  procedure.  DOT  is 
a  computer  program  created  by  VMA  Engineering  for  optimization.  DOT  uses  numerical 
optimization  to  change  the  selected  parameters  based  upon  the  search  gradients  to  either 
maximize  or  minimize  an  objective  function  nonlinear  in  the  design  variables.  The  main 
optimization  program  DOT  is  coupled  with  an  application  program,  which  determines  both  the 
constraints  and  the  objective  function,  by  writing  a  small  interface  FORTRAN  program  which 
controls  both  programs  [57],  An  outline  of  the  optimization  process  is  shown  in  figure  15. 
Prior  to  running  the  optimization,  realistic  material  data  on  the  piezoceramics  and  composite 
need  to  be  determined.  The  data  on  the  piezoceramic  is  provided  by  the  literature  from  the 
various  vendors  [43],  The  sensor  and  actuator  width,  wd,  is  0.01905  m  and  the  length  of  the 
actuator  is  0.0635  m  and  sensor  0.01905  m.  The  resistor  associated  with  an  ideal  operational 
amplifier  used  to  acquire  the  sensor  signal,  RF,  is  1E106  ohms.  The  elastic  modulus  of  both 
the  sensor  and  actuator  material,  cEu,  is  5.847  MPa.  The  actuator  location  in  the  thickness  or 
z-direction  is  set  equal  to  0.000635  m.  Finally,  the  piezoelectric  constant,  d31,  is  1.498E104 
m/V.  The  rest  of  the  constants  used  in  the  optimization  are  combinations  of  the  piezoelectric 
material  constants.  However,  the  composite  data  cannot  be  taken  directly  from  the  user 
material  data  sheet  because  the  material  ages  and  the  properties  are  distinct  for  each 
fabrication.  These  material  properties  must  be  determined  experimentally.  Tensile  testing  on 
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Figure  15  Outline  of  Optim  ization 


TEST  FOR  CONVEXITY 

f(pi) 

af(pl)  +  (1  -  a)f(p2)  < 

f(apl  +  (1  -  a)p2)  ri 

f(p2)  — 


apl  +  (1  -  a)p2  * 

P1 

P2 

0  <  a  <  1 


Figure  16  Test  for  Convexity 


three  sets  of  specimens  with  different  ply  angles  was  completed  to  determine  the  various 
composite  properties  required  for  this  research. 


The  procedure,  testing  method,  and  results  of  determining  the  composite  properties  are 
given  in  the  next  section  when  the  experimental  apparatus  is  discussed.  Only  the  results  of  the 
testing  are  presented  in  table  6. 


TABLE  6  EXPERIMENTALLY  DETERMINED  COMPOSITE  PROPERTIES 


Composite 

E, 

gi2 

V2| 

e2 

V12 

P 

Value 

156.45  GPa 

5.91  GPA 

.0107 

9.95  GPa 

.306 

1 

Prior  to  running  the  optimization  computer  program,  the  objective  function  was  tested 
for  convexity.  The  objective  function  and  the  constraint  equations  defined  as  the  functions 
f(p),  fl(p),  and  f2(p)  respectively. 

f(p)  =  t[Cr]Q[Cr]T  +  (1  -  C)  t Cb]Q[Cb]T  (154) 

flip)  =  cx  (155) 

£2(p)  =  C2  (156) 

In  this  research,  the  functions  fl  and  f2  pertain  to  each  actuator  and  sensor  pair.  Convexity 
was  tested  for  each  of  the  functions.  (See  figure  16.)  The  symbol  p  represents  a  complete  set 
of  optimization  parameters  such  as  the  ply  angle,  control  gains,  and  actuator  locations.  The 
symbol  a  is  an  arbitrary  multiplier  which  satisfies  0  <  a  <  1.  Some  randomly  generated 
optimization  parameters  within  the  appropriate  ranges  were  examined.  The  following  test  case 


variables  are  calculated  from  the  above  functions: 

TEST  =  af(pl)  +  (1  -  a)f(p2)  -  f(apl  +  (1  -  a)p2) 

TEST1  =  afl(pl)  +  (1  -  a)fl(p2)  -  fl(apl  +  (1  -  a)p2) 

TEST2  =  af2(pl)  +  (1  -  a)f2(p2)  -  f2(apl  +  (1  -  a)p2) 

With  this  setup,  convexity  was  tested  for.  If  any  of  the  TEST,  TEST1,  or  TEST2  variables 
were  negative  then  the  objective  function  or  the  constraint  equations  are  not  convex.  After 
thousands  of  runs,  the  objective  function  was  determined  to  be  not  convex,  because  there  were 
cases  in  which  the  TEST  function  was  negative.  This  means  that  a  local  minimum  may  not  be 
the  global  minimum  for  the  function.  However,  there  could  be  some  local  minima  which  can 
be  located  using  the  optimization  computer  program.  The  constraint  equations  are  considered 
to  be  convex,  because  in  every  case  both  the  TEST1  and  TEST2  function  were  positive.  This 
is  not  a  conclusive  test  for  convexity,  but  it  leads  one  to  believe  that  they  are  perhaps  convex. 
The  optimization  is  performed  in  order  to  determine  the  local  minima  and  try  to  find  the 
smallest  one  by  varying  the  optimization  parameters. 

In  order  to  use  DOT,  a  host  program  is  required  to  control  the  optimization  by 
defining  the  optimization  method  used,  number  of  design  variables  and  constraints,  and  by 
calling  the  main  optimization  program  as  a  subroutine.  This  shell  program  is  also  used  to 
define  the  programming  constants,  upper  and  lower  bounds,  and  initial  designs.  The  shell 
program  also  calls  a  subroutine  which  evaluates  the  objective  function  and  the  constraints  for 
the  optimization.  (See  equations  152-153.)  This  subroutine  takes  the  experimentally 
determined  composite  constants,  see  table  1,  and  the  sensor  and  actuator  constants  given  by 
equations  (80)  and  (77)  and  formulates  the  mass  and  stiffness  matrices  defined  by  equations 
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(104-105).  From  these  matrices,  the  eigenvalues  and  eigenvectors  are  determined  for  the  free 
motion  system.  Both  the  eigenvalues  and  eigenvectors  are  determined  by  using  subroutines 
from  the  IMSL  math  libraries.  Using  the  closed  loop  matrices,  the  Lyapunov  equation  (146)  is 
solved  and  used  to  calculate  the  performance  index,  equation  (152).  Finally,  the  constraint 
inequality,  equation  (153),  is  calculated.  All  of  these  programs  can  be  found  in  Appendix  B  of 
this  dissertation. 

In  this  research  the  optimization  method  used  by  DOT  is  Modified  Method  of  Feasible 
Directions.  This  method  is  reliable  and  uses  low  amounts  of  computer  memory.  The 
expressions  for  the  gradients  from  the  previous  theoretical  section  are  calculated  numerically 
by  the  finite  difference  method.  The  gradients  determine  which  search  direction  along  which 
the  optimization  program  is  to  proceed.  In  this  type  of  optimization,  the  process  moves  down 
to  the  minimum  in  the  direction  of  the  gradient  as  far  as  possible  until  a  constraint  is  reached. 
Once  a  constraint  is  reached,  a  new  search  direction  is  looked  for.  The  number  of  design 
variables  used  in  the  optimization  is  14:  10  angles  from  each  of  the  plies  in  the  composite 
layup,  2  from  the  control  gain  of  each  of  the  actuators,  and  2  from  the  x-location  of  each  of 
the  actuators.  Actually,  there  are  only  5  design  variables  for  the  ply  angles  due  to  the  forced 
symmetry.  For  each  of  the  ply  angles,  the  upper  bound  is  +90  degrees  and  the  lower  bound  is 
-90  degrees.  Anything  beyond  this  range,  the  material  properties  are  just  repeating  themselves. 
The  upper  bound  on  the  actuator  control  gains  is  40000  and  the  lower  bound  is  0.  This  upper 
bound  is  well  beyond  anything  we  could  produce  in  the  laboratory,  since  this  is  a  gain 
multiplier  which  amplifies  the  voltage  signal  from  the  sensor.  The  upper  and  lower  bounds  on 
the  actuator  x  axis  locations  are  0.5696m  and  0.0127m  which  correspond  to  the  beginning  and 
end  dimensions  of  the  plate.  The  complete  actuator  has  to  fit  on  the  plate  which  is  why  0  and 
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0.6096  were  not  used.  Finally,  £  was  set  equal  to  .5  which  corresponds  with  minimizing  the 
combined  bending  and  torsion  motion  of  the  plate.  The  parameters  that  were  held  fixed 
during  the  optimization  runs  were  the  piezoelectric  actuator  and  sensor  material  properties, 
composite  material  properties,  dimensions  of  the  plate,  location  of  the  sensors  and  actuators 
along  the  width  of  the  plate,  and  the  constant  The  constant  £  was  used  to  weigh  equally 
the  bending  tip  motion  with  the  torsional  tip  motion  of  the  plate.  The  main  objective  of  this 
research  is  to  stop  the  complete  motion,  both  bending  and  torsion,  of  the  plate  which  is  best 
accomplished  by  setting  C,  to  0.5.  (See  Table  7)  Later,  £  was  allowed  to  vary  from  0  to  1  by 
increments  of  .1  for  one  of  the  test  cases  to  show  the  affect  this  parameter  has  on  the 
optimization. 


TABLE  7  DEFINITION  OF  OPTIMIZATION  PARAMETERS 


Optimization  Parameter 

Held  Constant 

#  of  Variables  Optimized 

Ply  Angle,  9 

5 

Sensor/Actuator  X-Location 

2 

Control  Gains,  G 

2 

Actuator/Sensor  Properties 

X 

Composite  Properties 

X 

Plate  Dimensions,  a,b,t 

X 

Weighting  Parameter  C, 

X 

Sensor/Actuator  Y-Location 

X 

In  order  to  determine  if  the  minimum  objective  reached  during  the  optimization  routine 
is  the  lowest  local  minimum  or  possibly  the  global  minimum,  the  optimization  routine  was 
started  at  numerous  locations  for  the  design  variables.  The  starting  ply  was  varied  from  5  to 
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90  degrees  by  a  5  degree  differential  while  the  starting  gain  multiplier  is  varied  from  10  to 
.00001.  Both  starting  actuators  locations  were  varied  from  0.0254  to  0.55  meters  along  the 
length  of  the  plate.  The  results  of  the  optimization  parameter  study  can  be  found  in  Appendix 
C  of  this  dissertation.  Hundreds  of  computer  runs  were  made  to  determine  if  a  possible 
global  minimum  solution  could  be  found.  From  the  computer  runs,  the  value  of  the  objective 
function  reached  a  maximum  of  26875  which  occurred  with  a  starting  ply  angle  of  90  degrees, 
actuator  location  of  0.2032  and  0.4064  and  gain  multipliers  of  10.  The  final  results  of  this 
optimization  run  were  ply  angles  of  80.2  degrees,  actuator  location  of  0.188  and  0.577,  and 
actuator  gain  multipliers  of  8.86  and  0.12.  This  particular  result  from  the  optimization 
program  must  be  near  one  of  the  local  minima  since  the  optimized  parameters  are  located  near 
the  starting  parameters.  The  minimum  objective  function  of  approximately  60-70  was  reached 
consistently  by  about  50%  of  the  optimization  runs.  The  final  optimized  results  of  these 
optimization  runs  were  all  10  ply  angles  with  angles  of  about  20  degrees,  actuator  locations  of 
0.184  and  0.577  and  actuator  gains  of  8.23  and  0.142.  (See  Table  8  for  optimization  details) 


TABLE  8  SUMMARY  OF  "GLOBAL  MINIMUM"  SOLUTIONS  FROM  STARTING 

POSITIONS 


Parameters  to  be  Optimized 

Starting  Value 

Optimized  Value 

Ply  Angle  1 

5°  -  75° 

20° 

Ply  Angle  2 

5°  -  75° 

20° 

Ply  Angle  3 

5°  -  75° 

20° 

Ply  Angle  4 

5°  -  75° 

20° 

Ply  Angle  5 

5°  -  75° 

20° 

Ply  Angle  6 

5°  -  75° 

20° 

Ply  Angle  7 

5°  -  75° 

20° 

75 


Ply  Angle  8 

5°  -  75° 

O 

O 

<N 

Ply  Angle  9 

5°  -  75° 

20° 

Ply  Angle  10 

5°  -  75° 

20° 

Actuator  1  Location 

.2032  -  .4064  meters 

0.184  meters  from  base 

Actuator  2  Location 

.2032  -  .4064  meters 

0.577  meters  from  base 

Control  Gain  1 

.01  -  10.0 

8.23 

Control  Gain  2 

.01  -  10.0 

0.142 

Objective  Function 

60  -  70 

Based  upon  the  number  of  times  this  result  was  reached  no  matter  where  the  optimization  was 
started,  it  was  surmised  that  this  was  the  "global  optimal"  solution  for  this  problem. 

There  is  a  physical  significance  for  the  optimization  program  to  choose  those  particular 
parameters  as  the  global  solution.  Both  the  coupled  motion,  bending,  and  torsional  stiffness  of 
the  plate  are  being  used  to  stop  the  motion.  For  maximum  bending  stiffness,  the  optimum 
angle  is  0  degrees.  For  maximum  torsional  stiffness,  the  optimum  angle  is  45  degrees. 

Finally,  the  angle  which  maximizes  the  bending/torsional  coefficient,  D16,  is  1 8  degrees.  It  is 
reasonable  to  conclude  that  the  optimum  angle  would  occur  near  the  angle  which  maximizes 
this  coefficient.  It  turns  out  to  be  slightly  larger  due  to  the  fact  that  the  disturbance  is  a 
torsional  force  which  excites  predominately  only  the  torsion  and  bending/torsional  modes. 

In  order  to  verify  the  conclusions  concerning  the  optimized  ply  layup,  the  optimization 


program  is  rerun  with  different  performance  indices  as  C,  varies.  The  constant,  £,  is  varied 
from  0  to  1  by  increments  of  0.1.  The  results  are  presented  in  table  9. 


TABLE  9  EFFECTS  OF  £  ON  THE  PLY  ANGLE,  ACTUATOR  LOCATION  AND 

OBJECTIVE  FUNCTION 


G 

PLY  ANGLE 
(degrees) 

ACTUATOR.  - 
LOCATIONS 
(meters) 

OBJECTIVE  - 
FUNCTION 

0 

[_45/452/34/28]s 

.156/.577  from  base 

2.81 

0.1 

[262/9.8/21/16]s 

.186/.577  from  base 

21.3 

0.2 

[22/17/21/1 8/15]s 

.185/.577  from  base 

32.2 

0.3 

[203/19/16]s 

.185/.577  from  base 

42.2 

0.4 

[204/17]s 

.183/.577  from  base 

52.3 

0.5 

[205]s 

.184/.577  from  base 

61.7 

0.6 

[204/19]s 

.  1 837.577  from  base 

71.3 

0.7 

[1 9/2O3/1 8]s 

.  1 80/.577  from  base 

81.4 

0.8 

[204/17]s 

.169/.577  from  base 

92.2 

0.9 

[20,/17]s 

.161/.577  from  base 

101.3 

1.0 

[202/19/16]s 

.  1 497.577  from  base 

109.5 

From  table  9,  one  notices  that  the  objective  function  increases  linearly  as  £  changes  from  0  to 
1.  This  general  increase  in  the  objective  function,  with  increasing  £,  makes  physical  sense  due 
to  the  torsional  disturbance  input.  The  disturbance  excites  primarily  the  torsional  modes  of 
vibration  with  bending  coming  from  the  coupling.  One  would  expect  that  the  torsional  motion 
of  the  plate  would  be  larger  than  the  bending  motion  so  the  objective  function  would  increase 
as  the  dependance  upon  the  torsional  motion  is  increased.  One  should  note  that  for  the  C,  =  0 
case  which  is  controlling  pure  bending,  the  objective  function  should  approach  0  by  choosing  a 
balanced  layup  which  prevents  any  bending/torsional  coupling.  The  optimization  run  for  this 
case  is  tending  towards  a  balanced  layup  but  hits  a  local  minimum  or  simply  terminates 
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prematurely.  This  run  was  performed  for  only  one  starting  ply  and  actuator  location.  Also, 
notice  that  in  all  three  cases  the  actuator  location  does  not  vary  greatly. 

Finally,  it  is  important  to  determine  how  sensitive  the  optimized  objective  function  is 
to  small  changes  in  the  optimized  parameters.  If  the  ply  angle  changes  a  few  degrees,  or  if  the 
location  of  the  actuators  are  not  exactly  .184  and  .577  meters,  what  affect  does  this  have  on 
the  objective  function?  Using  slightly  different  parameters,  the  computer  program  which 
calculates  the  objective  function  when  supplied  the  optimized  parameters  was  run  again.  From 
these  results,  it  seems  clear  that  the  objective  function  is  not  sensitive  to  either  minor 
deviations  in  the  ply  angle  or  even  several  different  ply  angles  as  part  of  the  10.  In  fact,  in 
certain  cases  some  of  the  ply  angles  ended  up  at  40  degrees  and  the  objective  function  only 
changed  by  less  than  0.1%.  Also,  the  objective  is  not  sensitive  to  very  small  changes  in  the 
actuator  locations,  but  anything  larger  such  as  a  couple  of  centimeters  from  the  optimal 
locations  changes  the  objective  function  greatly.  If  the  actuator  location  varies  by  more  than 
0.0762  to  0.1016  meters,  then  the  objective  function  will  change  by  over  20%.  Also,  changing 
the  disturbance  force  does  not  change  some  of  the  optimization  parameters  such  as  the  actuator 
and  sensor  locations  and  ply  layup.  However,  the  change  in  the  disturbance  force  greatly 
effects  the  actuator  control  gains.  This  effect  is  almost  a  one-to-one  change.  If  the 
disturbance  force  is  doubled,  then  the  actuator  control  gains  are  reduced  by  one  half.  Also  if 
the  disturbance  force  is  halved,  then  the  actuator  control  gains  are  doubled.  This  relationship 
between  the  disturbance  force  and  the  actuator  control  gains  occurs  because  of  the  constraints 
placed  upon  the  total  amount  of  energy  supplied  to  the  actuators.  Because  of  physical 
limitations  of  the  actuator,  such  as  shorting  out  and  depoling,  the  total  voltage  supplied  to 
them  is  limited.  This  voltage  is  dependent  upon  both  the  strain  signal  of  the  sensor  and  the 
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actuator  control  gains.  In  fact  it  is  the  product  of  the  two.  Therefore  if  the  sensor  strain 
signal  increases,  then  the  actuator  control  gains  must  decrease  in  order  for  the  total  voltage  to 
remain  constant.  When  the  disturbance  force  is  increased,  the  strain  on  the  plate  is  also 
increased  which  causes  the  sensor  strain  signal  to  increase.  This  is  why  the  actuator  control 
gains  decrease  proportionally.  The  same  is  true  if  the  disturbance  force  is  decreased. 

The  next  chapter  contains  information  regarding  the  experimental  portion  of  this 
research.  A  plate  was  fabricated  and  tested  with  a  ply  angle  of  21  degrees.  The  results  of  this 
testing  were  then  compared  with  a  plate  fabricated  with  a  quasi-isotropic  layup. 
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V.  Experimental  Apparatus  and  Procedure 


The  information  presented  in  this  chapter  covers  the  experimental  apparatus  and 
procedure  for  the  fabrication,  the  testing  methodology,  the  material  characterization  results  of 
the  composite  test  specimens  and  the  testing  methodology  and  characterization  of  the  vibration 
testing.  The  procedures  used  to  characterize  the  material  follow  the  guidelines  provided  by 
ASTM.  The  standards  followed  for  determining  the  tensile  and  shear  properties  are  ASTM 
D3039-76  and  ASTM  D35 18-76.  The  size  of  the  specimen  and  the  testing  procedure  are  taken 
from  the  reference,  Experimental  Characterization  of  Advanced  Composite  Materials  by 
Carlsson  and  Pipes  [58],  The  specimens  were  fabricated  and  prepared  for  testing. at  the 
Composites  Structures  Laboratory  of  the  Applied  Composite  Branch,  Phillips  Laboratory 
Edwards  AFB,  CA.  The  material  used  was  Fiberite  Graphite/Epoxy  tape  or  specifically,  IM-7 
fiber  and  977-2  thermoset  epoxy  resin  on  12"  wide  preimpregnated  tape.  IM-7  fiber  is  a  high 
strength,  medium  modulus  graphite  fiber  that  has  been  used  extensively  for  space  applications. 
The  resin,  977-2,  is  a  dicyanate  ester  resin  which  is  space  qualified  and  provides  very  little 
microcracking  and  outgassing.  This  material  definitely  needed  to  be  characterized  because  it 
had  surpassed  its  shelf  life  by  about  2  years.  The  material  had  aged  significantly.  In  order  to 
completely  characterize  the  composite  material,  3  sets  of  samples  were  fabricated  with 
different  lay-ups  to  determine  all  of  the  material  properties.  The  samples  were  fabricated  with 
the  following  ply  angles:  [0]8  to  determine  E,  and  v12,  [90]16  to  determine  E2  and  v21,  and 
[±45]2s  to  determine  G12.  Also,  the  density  of  the  composite  material  was  determined  by 
measuring  the  volume  and  the  weight  of  some  of  the  test  specimens. 

The  fabrication  method  used  to  make  the  specimens  was  hand  layup.  This  is  the 
simplest  method  for  making  flat  panels  from  which  the  test  specimens  were  machined.  First, 
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the  composite  material  was  taken  out  of  a  freezer  and  allowed  to  thaw  while  the  mold  surface 
was  being  prepared.  This  preparation  involves  cleaning  the  mold  with  acetone  to  remove  any 
contaminates,  spraying  it  with  a  release  agent  such  as  McLube  or  Freecote,  and  then  attaching 
a  release  cloth  to  the  surface  using  teflon  tape  to  prevent  the  composite  from  adhering  to  the 
mold.  The  composite  material  used  was  .304  meters  wide  graphite/epoxy  unidirectional  tape 
which  means  the  reinforcing  fibers  all  run  in  only  one  direction.  The  composite  tape  was 
then  cut  into  square  and  rectangular  strips  depending  upon  which  layup  was  being  fabricated. 
These  strips  were  then  laid  on  top  of  one  another  in  the  proper  fiber  orientation  on  the  mold  to 
form  the  specimens.  A  combination  of  release  and  bleeder  cloth  is  then  attached  to  the  top  of 
the  mold  to  capture  the  excess  resin  which  flows  from  the  prepreg  material  out  of  the 
specimen  panel.  The  release  cloth  is  used  to  prevent  the  excess  resin  from  sticking  to  the 
formed  part.  A  vacuum  bag  is  then  formed  around  the  part  by  taping,  with  double-sided  tape, 
a  thick  plastic  sheet  to  the  mold.  (See  figure  17.)  Using  a  vacuum  pump,  the  air  is  removed 
from  the  bag  in  order  to  apply  uniform  pressure  on  the  composite  part  during  the  curing  cycle. 
The  part  is  now  ready  for  the  cure  cycle  which  completes  the  fabrication  of  the  composite 
material.  (See  figure  18)  The  cure  cycle  used  ramps  the  temperature  up  to  180  degrees  F.  At 
this  temperature  the  resin  becomes  less  viscous  and  starts  to  flow.  Once  this  temperature  is 
reached  80  psig  of  pressure  is  applied.  This  temperature  is  maintained  for  30  minutes  to  allow 
the  resin  to  heat  up  and  begin  to  flow  creating  a  uniform  cross-sectional  area  of  matrix 
material.  After  30  minutes  the  temperature  ramps  up  to  350  degrees  F  and  is  held  for  4  hours. 
The  resin  or  matrix  material  has  then  reached  its  TG,  or  gel  temperature,  and  begins  to  set. 

The  pressure  is  held  constant  during  this  time.  After  four  hours  both  the  temperature  and 
pressure  are  lowered  gradually  to  minimize  any  residual  thermostresses  which  could  be 
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STANDARD  CURE  CYCLE 
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TIME  (HOURS) 

Figure  18  Standard  Cure  Cycle  of  Graphic  Epoxy  Materials 
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induced  into  the  composite  material.  The  specimen  panels  are  now  complete  and  ready  to  be 
machined  to  test  specimens. 

Loading  tabs  are  attached  to  the  specimen  panels  in  order  for  the  material  testing 
machine  to  grip  the  individual  test  specimens.  These  loading  tabs  are  made  out  of  glass  epoxy 
and  are  attached  to  the  specimen  using  a  Hysol  9394  epoxy  resin.  The  resin  has  a  3-5  day 
cure  time  at  ambient  temperature  of  77  degrees  F  or  an  accelerated  cure  time  of  1  hour  at  150 
degrees  F.  The  loading  tabs  are  cut  into  strips  0.0381  meter  length  and  glued  to  the  specimen 
panels  prior  to  machining.  The  reason  the  loading  tabs  are  attached  prior  to  machining  is  to 
insure  that  the  specimen  and  tabs  are  aligned  with  one  another.  This  prevents  any  unwanted 
bending  force  from  entering  into  the  tensile  test.  After  the  resin  has  cured,  the  specimen 
panels  are  cut  into  their  proper  dimensions  depending  upon  the  ply  angle  orientation  of  the 
fiber  and  what  type  of  test  is  going  to  be  performed.  The  test  specimen  panels  were  cut 
according  to  the  specifications  stated  in  the  relevant  standards  for  shear  and  tension  testing, 
ASTM  D3039-76  and  D35 18-76.  Prior  to  testing,  electric  resistance  foil  strain  gages  are 
attached  both  in  the  longitudinal  and  transverse  directions  to  measure  the  strains  in  those 
directions.  These  gages  were  purchased  from  Measurements  Group,  Inc.  and  are  from  the 
CEA  -  Series.  They  are  self  compensating  strain  gages  with  a  fully  encapsulated  copper  grid 
and  exposed  copper-coated  integral  solder  tabs.  From  each  test,  different  material  constants 
are  derived.  The  specimens  were  loaded  into  a  MTS  material  testing  machine,  making  sure 
they  are  aligned  properly  to  prevent  bending  from  occurring.  The  specimens  are  loaded  at  a 
constant  rate  with  data  taken  continuously  until  failure.  Since  the  composite  material 
properties  are  directly  related  to  minute  changes  in  the  environment  during  the  fabrication  and 
curing  cycle,  numerous  tests  must  be  performed  and  then  averaged  to  obtain  the  material 
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constants.  Nine  specimens  of  each  ply  layup  were  tested  to  determine  a  statistical  average  of 
the  material  properties.  The  results  of  the  tests  for  0  degrees  are  shown  in  table  10. 


TABLE  10  RESULTS  FOR  THE  0  DEGREE  COMPOSITE  MATERIAL  TESTS 


Test  Number 

E, 

v,2 

1 

166.93  GPa 

Bad  Data 

2 

173.75  GPa 

Bad  Data 

3 

163.54  GPa 

.360564 

4 

174.71  GPa 

.38231 

5 

152.21  GPa 

.285173 

6 

159.14  GPa 

.307589 

7 

147.76  GPa 

.296064 

8 

127.49  GPa 

.263132 

9 

142.52  GPa 

.244482 

Average 

156.45  GPa 

.305616286 

The  results  of  the  tests  for  45  degrees  are  shown  in  table  1 1 . 
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TABLE  11  RESULTS  FOR  THE  45  DEGREE  COMPOSITE  MATERIAL  TESTS 


Test  Number 

Gu 

1 

5.4  GPa 

2 

5.65  GPa 

6.05  GPa 

4 

5.5  GPa 

5 

6.099  GPa 

6 

6.22  GPa 

7 

6.14  GPa 

8 

6.041  GPa 

Average 

5.91  GPa 

The  results  of  the  tests  for  90  degrees  are  shown  in  table  12. 


TABLE  12  RESULTS  FOR  THE  90  DEGREE  COMPOSITE  MATERIAL  TESTS 


Test  Number 

e2 

V2i 

1 

10.54  GPa 

.012422949 

2 

9.52  GPa 

.016004 

85 


3 

Bad  Data 

.018734 

4 

8.83  GPa 

.007410486 

5 

10  GPa 

Bad  Data 

6 

9.38  GPa 

0.004059426 

7 

10.61  GPa 

.008723356 

8 

10.2  GPa 

.009075485 

9 

9.72  GPa 

.0098146 

Average 

9.95  GPa 

.010780538 

The  density  was  determined  by  cutting  off  a  square  piece  of  composite  to  exact  dimension 
using  a  computer  controlled  milling  machine.  The  dimensions  are  checked  with  a  calibrated 
measuring  device.  Once  the  exact  volume  is  known,  the  specimen  is  weighed  on  a  precision 
scale.  The  density  is  then  calculated  from  these  experimental  measurements.  The  results  of 
these  characterization  tests  are  shown  in  table  6.  However,  as  one  can  see  there  is  a  lot  of 
variance  between  the  measured  values  of  v2l.  It  is  extremely  difficult  to  experimentally 
measure  v21,  because  the  strain  gage  slips  and  eventually  falls  off  prior  to  completing  the 
testing.  Prior  to  using  these  measured  numbers,  a  simple  check  was  made  using  the  equation: 


(157) 


It  was  determined  that  the  measured  value  for  v21,  0.010718,  was  inaccurate  and  the  number 
used  in  the  optimization  subroutine  was  calculated  to  be  0.01925.  A  similar  fabrication 
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procedure  described  above  to  make  the  test  specimens  was  used  to  fabricate  the  optimized 
plate  for  the  experiment. 

The  results  from  the  optimization  program  from  the  previous  chapter,  include  the 
optimal  ply  orientation  and  locations  of  the  actuators,  are  used  in  the  fabrication  of  the 
composite  plate.  This  structure  was  fabricated  at  the  Composite  Structures  Laboratory, 

Applied  Composite  Branch,  Phillips  Laboratory  at  Edwards  AFB,  CA.  The  plate  was  hand 
layed  up  on  an  aluminum  mandrel,  vacuum  bagged,  and  cured.  The  material  used  was  Fiberite 
Graphite/Epoxy  tape  or  specifically,  IM-7  fiber  and  977-2  Thermoset  Epoxy  Resin  on  0.304 
meter  wide  preimpregnated  tape.  Plates  having  a  quasi-isotropic  and  21  degree  layup  were 
fabricated.  Quasi-isotropic  means  that  the  properties  are  the  same  in  both  the  length  and  width 
direction.  This  is  accomplished  by  creating  a  layup  that  provides  both  the  stiffness  of  the 
fibers  and  the  elasticity  of  the  matrix  in  all  directions.  Plates  that  are  quasi-isotropic  usually 
have  a  layup  with  a  combination  of  0,  ±45,  and  90°  plies.  Fabricating  these  plates  with 
various  ply  orientations  is  extremely  difficult.  The  composite  material  is  similar  to  a  fabric 
with  sticky  glue  on  it  prior  to  its  being  cured.  Trying  to  cut  this  fabric  to  the  exact  geometry 
such  as  21  or  45  degrees  is  nearly  impossible.  Since  the  optimized  plies  varied  between  20 
and  22.05  degrees,  all  of  the  plies  of  the  optimized  plate  were  attempted  to  be  cut  at  21 
degrees.  The  difficulty  involved  in  cutting  causes  at  least  a  ±3  degree  error  in  the  ply  angle. 
Another  factor  which  influences  the  ply  angles  is  how  the  cut  composite  sheets  are  laid  on  top 
of  one  another.  When  the  plies  are  compacted  to  provide  adequate  bonding,  the  stretch  in  the 
matrix  directions  causes  the  ply  angle  to  change  somewhat.  These  errors  should  not  affect  the 
overall  results  of  this  research  since  it  was  shown  in  the  last  chapter  that  the  optimized 
objective  function  is  not  extremely  sensitive  to  changes  in  the  ply  orientation.  The  plates  were 
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cured  using  a  standard  cure  cycle  at  176.67  degrees  C  at  551.6  kPa  of  pressure,  the  same  as 
the  test  coupons.  After  the  plates  were  fabricated,  the  piezoceramic  actuators  and  sensors  were 
attached  to  the  plate  at  the  optimized  locations  as  determined  by  the  optimization  computer 
program.  The  sensors  and  actuators  were  attached  using  room  temperature  cure  epoxy  resin. 
The  sensors  are  attached  on  one  side  of  plate  and  the  actuators  are  attached  in  the  exact 
location  on  the  opposite  side  of  the  plate.  (See  figure  19.)  Once  finished,  this  completed  plate 
is  then  attached  to  the  control  electronics. 

The  structural  assembly  described  above  is  attached  to  the  control  electronics  described 
in  the  theoretical  section  which  contains  among  other  things,  a  strain  rate  feedback  circuit. 

This  circuit  is  shown  in  figure  20.  For  the  testing,  two  circuits  were  built  out  of  breadboard 
electronics.  The  responses  of  both  circuits  targeted  the  torsional  mode  of  the  plate.  This 
means  that  the  circuits  provided  the  most  damping  at  the  torsional  theoretical  resonance 
frequency  of  the  plate.  This  response  is  shown  in  figures  21,  22.  The  output  of  the 
piezoceramic  sensors  is  conditioned  by  the  feedback  electronics  which  then  supply  a  high 
voltage  control  signal  to  an  actuator.  One  sensor  will  control  only  one  actuator  in  this 
research.  The  feedback  electronics  consist  of  a  preamplifer  which  converts  the  sensor  signal  to 
a  voltage,  a  bandpass  filter  with  an  adjustable  center  frequency  and  bandwidth  which  changes 
the  signals  phase,  and  a  high  power  output  amplifier  which  amplifies  the  voltage  signal  to 
power  the  actuators.  The  center  frequency  of  the  electronic  circuit  is  centered  around  the 
torsional  and  bending/torsional  modes  of  the  plates.  The  torsional  vibration  is  provided  by 
using  a  piezoceramic  sensor  and  actuator  pair  mounted  halfway  along  the  length  of  the  plate  at 
a  45  degree  to  the  longitudal  axis.  A  pseudo-random  signal  of  500  mVrms  is  sent  to  the 
actuators  as  the  disturbance  source. 


88 


Location  of  Actuator  and  Sensor  on  Plate 


Figure  19  Location  of  Sensor  and  Actuator  on  Plate 


Strain  Rate  Feedback  Circuit 


Figure  20  Strain  Rate  Feedback  Circuit 
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RESPONSE  -  STRAIN  RATE  FEEDBACK  CIRCUIT  1 


Harken  X  Ref  A:  O  X  Ref  S:  O 

Y  Ref  A:  o  Y  Ref  B:  o 

Date:  OS-12-S4  Time:  OB:  22:  OO  AM 


Figure  21  Response  -  Strain  Rate  Feedback  Circuit  1 


RESPONSE  -  STRAIN  RATE  FEEDBACK  CIRCUIT  2 


Source  Level:  500  mVrme 

[FFT] 

Date:  00-12-94  Time:  08:  37:  00  AH 


A  schematic  of  the  control  electronics  is  shown  in  figure  23. 

The  main  objective  of  the  experimental  portion  is  to  characterize  the  damping  of  the 
plate  and  its  performance  when  subjected  to  random  noise  disturbances.  The  structural 
response  can  be  measured  by  exciting  the  structure  with  a  piezoceramic  actuator  and 
measuring  the  response  with  a  piezoceramic  sensor.  The  complex  modal  parameters  and  mode 
shapes  of  the  structure  will  be  calculated  using  transfer  function  analysis  techniques  as 
determined  by  the  frequency  response  functions  (FRF).  Initially,  the  plates  are  characterized 
with  the  control  system  turned  off.  The  random  noise  signal  is  sent  to  the  disturbance 
piezoceramic  actuator  and  the  structural  response  is  detected  by  the  collocated  sensor.  The 
response  is  the  real  time  average  of  100  signals.  Once  the  plates  are  characterized,  then  the 
above  procedure  is  repeated  with  the  control  system  turned  on.  The  closed  loop  response  of 
the  21  degree  plate  is  measured.  The  magnitude  of  the  natural  frequencies  are  compared  with 
the  theoretical  predictions.  A  HP  35665A  Hewlet  Packard  Dynamic  Signal  Analyzer  is  used 
to  create  the  random  noise  signal  and  to  calculate  the  (FRF)'s  of  the  plates.  The  signal 
analyzer  is  also  used  to  characterize  the  response  of  each  of  the  electronic  circuits  used  to 
control  the  plate. 
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ACTUATOR  ELECTRICAL  PATHWAYS 
SENSOR  ELECTRICAL  PATHWAYS 


Figure  23  Schematic  of  Control  Electronics 


VI.  Comparisons  Between  Baseline,  Modified  Plate  and  Theory 


In  this  chapter,  the  21  degree  optimized  plate  is  compared  to  a  quasi-isotropic  baseline 
plate  as  well  as  the  theoretical  calculations  in  terms  of  performance.  -When  comparing  the 
optimized  plate  to  the  theoretical  calculations,  one  should  note  that  the  actuator  locations  for 
the  experimental  plate  are  not  at  the  optimized  locations.  They  were  placed  at  their  current 
location  of  0.531m  based  upon  theoretical  calculations  that  were  in  error.  That  error  has  been 
corrected  producing  the  new  actuator  locations,  however  a  new  experimental  plate  was  not 
fabricated.  One  should  also  note  that  the  optimal  ply  angles  did  not  change  due  to  the  error  in 
the  theory.  The  theoretical  calculations  and  experiments  are  also  compared  with  the  work 
done  by  Hwang,  Hwang,  and  Chul  [25]  [26],  Their  work  is  very  similar  to  the  work 
accomplished  here  with  the  exception  of  their  modelling  method.  As  stated  previously,  the 
results  found  in  Hwang,  Hwang,  and  Chul's  report  are  very  close  to  the  results  found  in  this 
dissertation  even  though  the  modelling  methodology  is  completely  different:  finite  element 
method  versus  classical  Rayleigh-Ritz  modelling  and  sequential  versus  simultaneous 
optimization.  First,  the  work  accomplished  in  the  above  paper  is  compared  to  theoretical 
calculations  in  this  dissertation. 

In  the  papers  by  Hwang,  Hwang,  and  Chul  [25]  [26],  numerical  calculations  are 
performed  on  an  eight-ply  laminated  composite  plate  with  two  of  the  plies  variable.  The 
stacking  sequence  is  [9|/92/0790°]s  where  9,  and  92  are  the  ply  angle  design  variables.  The 
optimization  is  first  conducted  for  an  orthotropic  plate  with  the  design  angles  set  equal  to  0 
degrees  to  determine  the  optimal  size  and  location  of  the  piezoelectric  sensors  and  actuators. 
Initially,  the  size  of  the  sensors  and  actuators  are  unconstrained  while  the  locations  are  not 
allowed  to  vary.  In  the  second  case  the  actuator  sizes  are  constrained  while  the  optimal 
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locations  are  obtained.  The  optimization  conducted  next  is  to  determine  the  optimal  ply  layer 
angles  of  the  laminated  plate  which  will  maximize  the  control  performance.  The  optimization 
is  conducted  with  and  without  dynamic  stability  constraints.  The  optimization  is  performed  for 
only  a  limited  number  of  cases.  The  results  of  the  optimized  actuator  locations  are  shown  in 
figure  24.  As  the  reader  can  see  they  compare  favorably  with  what  this  dissertation  predicts. 
The  theory  in  this  dissertation  places  the  actuators  at  the  optimum  location  of  0.184  and  0.577 
meters  from  the  base  out  of  a  plate  length  of  0.6096  meters  or  30%  and  95%  down  the  length 
of  the  plate.  In  Hwang's  paper,  the  optimum  location  of  the  actuators  are  ,41m  from  the  base 
out  of  a  plate  length  of  ,5m  or  82%  down  the  length  of  the  plate.  Thus  there  is  only  a  13% 
difference  between  one  actuator  location  in  these  two  theories  which  were  developed 
differently  and  independently.  The  results  of  the  optimized  ply  angles  also  compare  very 
favorably  to  the  theoretical  predictions  in  this  dissertation.  These  results  are  shown  in  table 
13.  The  results  from  this  dissertation  predicts  the  ply  angles  to  be  20-21  degrees  while  the 
results  from  HHC's  paper  vary  from  24.98  -  29.04  degrees.  These  ply  angles  differ  from  this 
dissertation  by  20  -  31%  depending  upon  which  ply  angle  you  choose.  This  may  seem  high, 
but  remember  that  only  two  of  the  ply  angles  are  allowed  to  vary  in  Hwang's  paper  and  the 
others  are  fixed.  Since  the  fixed  ply  angles  offer  very  little  torsional  stiffness,  the  optimized 
ply  angles  are  expected  to  be  higher  or  closer  to  the  45  degree  maximum  torsional  stiffness. 
Table  14  shows  the  final  optimized  results  of  the  two  theories. 
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TABLE  13  OPTIMAL  PLY  ANGLE  FROM  [24]  [25]  [26] 


Objective  Function 

Case  No. 

Design  Variable  01 , 

,  Design  Variable  02 

C3  for  IB 

1 

30.00/26.19 

30.00/28.75 

C3  for  IB 

2 

30.00/25.10 

-30.00/26.10 

C3  for  IT 

1 

30.00/26.12 

30.00/29.04 

C3  for  IT 

2 

30.00/25.94 

-30.00/24.98 

TABLE  14  COMPARISON  BETWEEN  REFERENCED  PAPER  AND  DISSERTATION 

PARAMETERS 


Optimized  Parameter 

Referenced  Paper  Theory 

Dissertation  Theory 

Actuator  Location  #1 

,41m/.5m  82% 

0.184m/0. 6096m  30% 

Actuator  Location  #2 

,41m/.5m  82% 

0.577m/0. 6096m  94% 

Ply  Angle  #1 

24.98  -  29.04  degrees 

20-21  degrees 

Ply  Angle  #2 

24.98  -  29.04  degrees 

20-21  degrees 

The  results  of  the  experiment  will  now  be  presented  and  compared  with  the  theoretical 
calculations.  The  theoretical  modes  of  vibration  are  shown  in  figures  25  -  28.  Each  figure 
displays  a  different  mode  of  vibration:  figure  25  -  1  st  bending,  figure  26  -  2nd  bending,  figure 
27  -  1st  torsion,  figure  28  -  2nd  torsion.  The  frequencies  for  the  theoretical,  both  4  and  9 


mode  models,  and  experimental  calculations  of  the  21  degree  plate  are  shown  in  table  15. 


TABLE  15  FREQUENCY  COMPARISON  OF  THEORETICAL  CALCULATIONS  AND 
EXPERIMENTAL  MEASUREMENTS  OF  THE  21  DEGREE  PLATE 


Frequency 

4  Mode 

9  Mode 

Exp 

4  Mode  vs  Exp 

9  Mode  vs  Exp 

1st  Bending 

1.224 

1.15 

1 

22% 

15% 

2nd  Bending 

8.62 

7.00 

6 

44% 

17% 

1st  Torsion 

16.35 

14.66 

18 

-9.2% 

-19% 

2nd  Torsion 

64.54 

53.10 

54 

20% 

-1.7% 

The  percent  difference  in  the  frequencies  between  both  the  theoretical  calculations  and  the 
experimental  measurements  are  not  that  significant  when  modelling  composite  structures 
especially  with  the  nine  mode  theoretical  model.  As  one  can  see,  in  all  but  one  case  the 
Rayleigh-Ritz  approximate  solution  provides  an  upper  bound  to  the  experimental  frequencies. 
One  reason  that  theoretically  determined  and  experimentally  measured  frequencies  varied  is  the 
added  weight  and  stiffness  associated  with  the  piezoceramics  were  not  accounted  for  in  the 
analysis.  The  plate  itself  weighed  only  149.9  g,  each  of  the  three  actuator  elements  weighed 
4.8  g,  and  each  of  the  three  sensor  elements  weighed  2.4  g  bring  the  total  weight  of  the 
structure  to  171.5  g.  This  is  an  increase  in  weight  of  14.4%.  The  actuators  and  sensors  also 
provide  extra  stiffness  in  certain  modes  of  the  structure  such  as  the  torsional  mode. 

As  far  as  how  the  optimized  plate  performed  experimentally,  both  the  open  and  closed 
loop  response  to  the  21  degree  optimized  plate  showed  a  significant  improvement  in  stopping 
the  tip  motion  over  the  open  loop  baseline  isotropic  plate.  The  performance  of  the  open  loop 
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response  of  the  isotropic  plate  is  shown  in  figure  29.  The  performance  comparison  between 
the  open  and  closed  loop  response  of  the  optimized  plates  are  shown  in  figure  30.  Figures  29 
and  30  show  that  both  the  open  and  closed  loop  response  of  the  optimized  plate  perform 
significantly  better  than  the  open  loop  response  of  the  isotropic  plate.  The  open  loop  response 
of  the  optimized  plate  was  16.6%  lower  than  the  open  loop  response  of  the  isotropic  plate. 

The  closed  loop  response  performed  even  better.  The  closed  loop  response  of  the  optimized 
plate  was  29.8%  lower  than  the  open  loop  response  of  the  isotropic  plate. 

In  the  case  of  the  optimized  plate,  the  performance  improved  with  the  addition  of  the 
control  system.  The  magnitude  of  the  harmonic  peak  for  the  1st  torsional  mode  decreased  by 
11.3%  and  for  the  2nd  torsional  mode  by  6.7%.  Both  the  isotropic  and  theoptimized  plate 
tests  were  repeated  10  times  with  excellent  repeatable  results.  One  comparison  that  was 
unable  to  be  made  due  to  the  limitations  with  the  equipment  that  was  used,  was  comparing  the 
damping  factors  of  the  theoretical  predictions  with  the  experimental  measurements.  The 
damping  factors  are  measurements  of  the  inherent  structural  damping  or  added  damping  due  to 
the  active  control  system.  Next,  the  changes  in  the  damping  factors  for  the  different  ply  are 
presented.  Note,  the  actuator  gains  and  locations  are  the  same  for  both  plates  in  this 
comparison.  These  are  shown  in  table  16.  These  were  determined  from  the  theoretical 
analysis  from  the  closed  loop  plant  matrix. 
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TABLE  16  THEORETICAL  DAMPING  FACTORS  FOR  ISOTROPIC  AND  21  DEGREE 


PLATE 


Mode  of  Vibration 

Isotropic  Plate 

21  degree  Optimized  Plate 

1st  Bending 

.00039 

.000248 

2nd  Bending 

.0945 

.025 

1st  Torsion 

.0129 

.0685 

2nd  Torsion 

.0000009781 

.000084936 

The  damping  factor  is  greater  in  the  isotropic  plate  than  in  the  optimized  plate  for  the  bending 
modes  of  vibration.  However,  as  one  would  expect,  the  damping  factors  are  greater  in  the 
torsion  and  combined  bending/torsion  modes  of  vibration.  In  fact,  for  the  first  torsional  mode 
the  difference  is  a  factor  of  5  and  for  the  second  torsional  mode  the  difference  is  two  orders  of 
magnitude.  This  implies  that  the  isotropic  plate  damps  out  the  bending  modes  of  vibration 
better  than  the  optimized  plate.  However,  the  optimized  plate  damps  out  both  the  1st 
torsional  and  2nd  torsional  modes  of  vibration  better  than  the  isotropic  plate.  Since  the 
optimized  plate  was  designed  specifically  to  stop  the  motion  of  the  plate  when  subjected  to  a 
torsional  disturbance,  these  results  are  expected. 
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VII.  Summary 


The  results  from  this  research  are  a  design  methodology  and  a  plate  which  is  optimized  with 
respect  to  both  structural  and  control  parameters.  The  design  methodology  demonstrates  how  to 
achieve  simultaneous  structural  and  control  optimization  using  the  basic  elastic  and  dynamic  equations. 
The  ply  layup  and  location  of  the  piezoelectric  actuators  that  provide  bending  torque  only  are 
determined  from  the  optimization  routine.  The  ply  orientation  of  the  plate  is  optimized  to  enhance  the 
active  control  elements.  The  active  control  system  is  optimized  with  respect  to  actuator  and  sensor 
location  and  the  amount  of  voltage  being  applied  to  each  of  the  actuators.  Finally,  the  optimized  ply 
orientation  produces  coupling  between  the  bending  and  torsion  modes  of  vibration.  With  this 
"tailored"  coupling  not  only  will  the  piezoelectric  actuators  damp  bending,  but  they  will  also  damp 
the  rotation  motion. 

The  theoretical  calculations  provided  the  natural  frequencies,  mode  shapes  and  damping  factors 
for  the  plates  which  corresponded  very  well  with  the  measured  experimental  values.  The  percent 
difference  of  the  natural  frequencies  varied  by  only  9%  to  44%  which  is  well  within  experimental 
error  when  modelling  composite  structures.  The  theory  also  predicted  that  the  optimized  plate  would 
achieve  higher  damping  in  the  torsional  mode  than  the  quasi-isotropic  plate.  The  factor  of  5  increase 
in  torsional  damping  is  a  significant  increase  in  damping  considering  the  fact  that  the  method  used 
secondary  effects.  The  torsion  modes  are  not  being  damped  directly  as  in  the  1st  or  2nd  bending 
modes.  Recall  that  the  piezoceramic  actuators  only  produce  a  bending  moment  which  easily  reacts 
against  the  bending  motion  of  the  plate.  This  direct  stopping  motion  produces  large  damping  in  the 
1st  or  2nd  bending  modes  of  the  plate.  Only  the  bending/torsional  coupling  motion  of  the  plate 
caused  by  the  optimized  ply  layup  and  actuator  location  produces  any  increase  in  damping  for  a 
torsional  disturbance  force.  This  is  verified  by  the  figures  in  chapter  VI,  which  show  both  the  open 
and  closed  loop  response  for  the  baseline  plate  and  the  optimized  plate. 
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The  present  theoretical  calculations  compare  favorably  with  independent  papers  written  by 
Hwang,  Hwang,  and  Chul  [25]  [26].  The  concept  of  their  papers  and  referenced  PhD  dissertation  by 
W.  S.  Hwang  in  1993,  is  almost  identical  to  what  is  proposed  in  this  dissertation.  However,  it  was 
shown  that  in  both  chapters  II  and  VI,  their  theory  and  approach  are  different.  Their  theory  is  based 
upon  finite  element  modelling  and  their  optimization  subroutine  is  performed  sequentially  whereas  this 
dissertation  models  the  plate  by  using  energy  method  models  and  simultaneous  optimization.  An  even 
greater  difference  between  the  two  theories  is  in  the  geometry  of  the  plate  and  materials  used.  Also, 
experiments  were  conducted  in  this  dissertation  to  back  up  the  theory,  while  the  referenced  papers 
discuss  only  theoretical  work.  Even  though  the  differences  in  the  theory  are  significant,  what  is  even 
more  outstanding  is  the  fact  that  the  results  compare  so  well.  One  of  the  optimized  actuator  locations 
differ  by  only  13%  and  the  optimized  ply  angles  differ  by  20%  to  31%.  Based  upon  the  results  of  this 
dissertation  and  the  other  referenced  work,  the  concept  of  damping  out  torsional  vibrations  using 
bending  actuators  enhanced  by  the  coupled  bending  torsional  motion  of  a  plate  is  achievable  and 
worthwhile. 

Further  refinement  of  this  theory  with  the  addition  of  a  combined  bending/torsional  disturbance 
force  and  taking  into  account  the  viscoelastic  damping  effects  of  the  matrix  material  would  be  a  good 
topic  for  further  research.  Substituting  a  combined  bending/torsional  force  as  the  disturbance  instead 
of  a  torsional  force  would  more  than  likely  change  the  ply  orientation  and  at  least  one  of  the  actuator 
locations.  These  would  change  in  order  to  provide  directly  damping  to  the  induced  bending  modes  of 
vibration  of  the  plate.  Also,  taking  into  account  the  inherent  damping  effects  of  the  composite 
material  would  show  conclusively  the  effects  of  the  active  control  system.  One  other  minor 
modification  which  could  be  made  to  refine  this  theory  would  be  to  allow  some  of  the  fixed  structural 
parameters  to  vary  such  as  the  y  -  direction  location  of  the  actuator/sensor  pair.  If  this  theory  and 
design  approach  was  used  in  an  operational  satellite  system,  the  plate  length  and  width  could  also  be 
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allowed  to  vary  to  determine  the  complete  optimized  structure  for  a  given  disturbance  profile. 

As  stated  in  the  introduction,  as  satellites  require  more  power,  the  size  of  the  solar  array 
increases  tremendously  which  causes  the  solar  arrays  to  be  prone  to  vibrations  due  to  their  large  size 
and  light  weight.  A  design  methodology  and  structure  fully  optimized  with  respect  to  its  size,  shape, 
and  structural  makeup  which  can  damp  out  both  bending  and  torsional  vibration  would  be  greatly 
beneficial  to  the  military  and  commercial  space  community.  This  research  centered  on  developing  and 
testing  this  design  methodology  by  fabricating  an  optimized  plate  which  provides  positive  proof  of  the 
theoretical  calculations. 


105 


Bibliography 


1.  Agrawal,  B.  N..  Design  of  Geosynchronous  Spacecraft.  Prentice-Hall,  Inc.,  Englewood 
Cliffs,  NJ,  1986. 

2.  Forrester,  P.  G.,  "Effects  of  Aeroelastic  Tailoring  on  Anisotropic  Composite  Material 
Beam  Models  of  Helicopter  Blades,"  Master's  of  Science  Thesis  University  of  Virginia, 
Richmond,  VA,  May  1989. 

3.  Jones,  R.  M.,  Mechanics  of  Composite  Materials.  Hemisphere  Publishing  Corp.,  New 
York,  1975. 

4.  Murray,  Y.  D.,  "Fibrous  Composite  Laminates  -  An  Introduction  to  the  Theory  and 
Verification  of  the  DYNA3D  Implementation,"  Technical  Report  #APTEK,  Aug  1988. 

5.  Cook,  R.  D.,  Malkus,  D.  S.,  and  M.  E.  Plesha,  Concepts  and  Applications  of  Finite 
Element  Analysis.  Wiley  Publishing,  New  York,  1989. 

6.  Weisshaar,  T.  A.,  "Aeroelastic  Tailoring  -  Creative  Uses  of  Unusual  Materials," 
AIAA-87-0976  AIAA  Paper,  1987. 

7.  Shirk,  M.  H.,  Hertz,  T.  A.,  and  T.  A.  Weisshaar,  "Aeroelastic  Tailoring  -  Theory, 
Practice,  Promise,"  Journal  of  Aircraft.  Vol  23,  No.  1,  pp.  6-18,  Jan  1986. 

8.  Crawley,  E.  F.,  and  K.  B.  Lazarus,  "Induced  Strain  Actuation  of  Isotropic  and 
Anisotropic  Plates,"  AIAA  Journal.  Vol.  29,  No.  6.,  pp.  944-951,  June  1991. 

9.  Fukunaga,  H.  and  H.  Sekine,  "A  Laminate  Design  for  Elastic  Properties  of  Symmetric 
Laminates  with  Extension-Shear  or  Bending-Twisting  Coupling,"  Journal  of  Composite 
Material.  Vol.  28,  No.  8,  pp.708-731,  1994. 

10.  Collar,  A.  Ri,  "The  First  Fifty  Years  of  Aeroelasticity,"  Aerospace.  Vol  5,  No.  2,  pp. 
111-125,  May  1978. 

11.  -  Ashley,  H.,  "The  Constructive  Uses  of  Aeroelasticity,"  AIAA-80-0877  AIAA  Paper, 
1982. 

12.  Garfinkle,  M.,  "Twisting  Smartly  in  the  Wind,"  Aerospace  America.  Vol.  32,  No.  7, 
pp.  18-20,  1994. 

13.  Lake,  R.  C.,  and  M.  W.  Nixon,  "A  Preliminary  Investigation  of  Finite-Element 
Modeling  for  Composite  Rotor  Blades.  Aerostructures  Directorate."  Proceeding  of  the  3rd 
Workshop  on  Dynamics  of  Aeroelastic  Structures.  NASA  Langley  Research  Center,  Virginia, 
Vol.  1,  pp.  354-365,  1990. 


106 


14.  Weisshaar,  T.  A.  and  S.  M.  Ehlers,  "Adaptive  Aeroelastic  Composite  Wings-Control 
and  Optimization  Issues,"  Journal  of  Composite  Engineering,  Vol.  2,  No.  5-7,  pp.  457-476, 
1992. 


15.  Hodges,  D.  H.,  Atilgan,  A.  R.,  Fulton,  M.  V.  and  L.  W.  Rehfield,  "Free- Vibration 
Analysis  of  Composite  Beams."  Proceedings  of  American  Helicopter  Society  National 
Specialists  Conference.  Vol.  1,  pp.  36-47,  Arlington,  TX,  1989. 

16.  Hwang,  S.  J.,  and  R.  F.  Gibson,  "The  Influence  of  Vibration  Coupling  on  Damping  of 
Laminated  Composites,"  AIAA  Journal.  Vol.  29,  No.  10,  pp.  1678-1685,  Oct  1991. 

17.  Barrett,  D.  J.,  "A  Design  For  Improving  the  Structural  Damping  Properties  of  Axial 
Members,"  Proceedings  of  Damping  1989.  WRDC-TR-89-31 16.  Vol.  Ill,  pp.  HCB-l-HCB-20, 
8-10  February  1989. 

18.  Rotz,  C.  A.  and  D.  J.  Barrett,  "Cocured  Damping  Layers  in  Composite  Structures," 
Proceedings  23rd  International  SAMPE  Technical  Conference.  Vol.  23,  pp.  352-362,  Anaheim, 
CA,  Oct  1991. 

19.  Belknap,  F.  M.  and  J.  B.  Kosmatka,  "Vibration  Suppression  of  Thin-Walled  Composite 
Tubes  using  Embedded  Viscoelastic  Layers,"  Proceedings  -  Damping  1991.  pp.  HAC-1  - 
HAC-16,  San  Diego,  CA,  Feb  1991. 

20.  Bronowicki,  A.  J.  and  H.  P.  Diaz,  "Analysis,  Optimization,  Fabrication,  and  Test  of 
Composite  Shells  with  Embedded  Viscoelastic  Layers,"  Proceedings  of  Damping  1989. 
WRDC-TR-89-31 16,  Vol.  Ill,  pp.  GCA- 1 -GCA-21,  8-10  February  1989. 

21.  Olcott,  D.  D.,  Improved  Damping  in  Composite  Structures  Through  Stress  Coupling. 
Co-cured  Damping  Lavers,  and  Segmented  Stiffness  Lavers.  Ph.D.  Dissertation,  Brigham 
Young  University,  Salt  Lake  City,  UT,  1992. 

22.  Olcott,  D.  D.,  Rotz,  C.  A.,  and  D.  J.  Barrett,  "Improved  Damping  in  Composite  Tubes 
Through  Stress  Coupling  and  Co-cured  Damping  Layers,"  Proceedings  23rd  International 
SAMPE  Technical  Conference.  Vol.  23,  pp.  373-387,  Anaheim,  CA,  Oct  1991. 

23.  *  Olcott,  D.  D.,  Rotz,  C.  A.,  and  P.  F.  Eastman,  "Improved  Vibration  Damping  in 
Composite  Structures  Using  'Zig-Zag'  Fibers  and  Embedded  Viscoelastic  Damping  Layers," 
Proceedings  38th  International  SAMPE  Symposium.  Vol.  1,  pp.  1357-1370,  Anaheim,  CA, 

May  1993. 

24.  Hwang,  W.,  Park,  H.  C.,  and  W.  Hwang,  "Tailoring  of  a  Laminated  Composite  Plate 
Under  Structural  Control,"  Proceedings  of  ICCM/9:  Composite  Modelling  and  Processing 
Science.  Madrid,  Spain,  Vol.  3,  pp.  592-598,  12-16  July  1993. 

25.  Hwang,  W.;  Park,  H.  C.;  and  W.  Hwang,  "Vibration  Control  of  Laminated  Plate  with 
Piezoelectric  Sensor/Actuator:  Finite  Element  Formulation  and  Modal  Analysis,"  Journal  of 


107 


Intelligent  Material  Systems  and  Structures.  Vol  4,  pp.  317-329,  July  1993. 


26.  Hwang,  W.;  Hwang,  W.,  and  H.  C.  Park,  "Integration  of  Composite  Structural  Design 
with  the  Intelligent  System  Concept,"  Proceedings  of  the  34th  AIAA/ASME/ASCE/AHS/ASC 
Structures.  Structural  Dynamics,  and  Materials  Conference,  part  6  pp.  3534-3539,  1993. 

27.  Shen,  I.  Y.,  "Bending  and  Torsional  Vibration  Control  of  Composite  Beams  Through 
Intelligent  Constrained-Layer  Damping  Treatment,"  Proceedings  of  1995  Smart  Structures  and 
Materials  Conference.  San  Diego,  CA,  Vol.  2445,  pp.  110-122,  Mar  1-2  1995. 

28.  Napolitano,  K.  L.  and  J.  B.  Kosmatka,  "Extension-Twist  Coupled  Damped  Composite 
Strut:  Analytical  and  Experimental  Evaluation,"  Proceedings  of  1995  Smart  Structures  and 
Materials  Conference.  San  Diego,  CA,  Vol.  2445,  pp.  362-373,  Mar  1-2  1995. 

29.  Humphreys,  E.  A.,  and  B.  W.  Rosen,  "Passive  Damping  Characteristics  of  Satellite 
Structural  Components",  WRDC  Technical  Report,  WRDC-TR-89-4048,  July  1989. 

30.  Adali,  S.,  and  V.  Verijenko,  "Minimum  Weight  Design  of  Symmetrical  Angle-Ply 
Laminates  Under  Multiple  Uncertain  Loads,"  Journal  of  Structural  Optimization.  Vol.  9,  No.  2, 
pp.  89-95,  Apr  1995. 

31.  Walker,  M.,  Adali,  S.,  and  V.  Verijenko,  "Optimization  of  Symmetric  Laminates  for 
Maximum  Buckling  Load  Including  the  Effects  of  Bending-Twisting  Coupling,"  Journal  of 
Computers  and  Structures.  Vol.  58,  No.  2,  pp.  313-319,  Jan  1996. 

32.  Dracopoulos,  T.  N.,  and  H.  Oz,  "Integrated  Aeroelastic  Control  Optimization  of 
Laminated  Composite  Lifting  Surfaces,"  Journal  of  Aircraft.  Vol.  29,  No.  2,  pp.  280-288, 
March-April  1992. 

33.  Szilard,  R.,  Theory  and  Analysis  of  Plates.  Prentice-Hall,  Inc.,  Englewood  Cliffs,  NJ, 
1974. 

34.  Fung,  Y.  C.,  Foundations  of  Solid  Mechanics.  Prentice-Hall  Inc.,  Englewood  Cliffs, 

NJ,  1965. 

35.  -  Meirovitch,  L.,  Analytical  Methods  in  Vibration.  Macmillan  Publishing  Co.  Inc.,  New 
York,  NY,  1967. 

36.  Clark,  S.  K.,  Dynamics  of  Continuous  Elements.  Prentice  Hall  Inc.,  Englewood  Cliffs, 
NJ,  1972. 

37.  Meirovitch,  L.,  Methods  of  Analytical  Dynamics.  McGraw-Hill  Publishing  Company, 
New  York,  NY,  1970. 

38.  Crawley,  E.  F.  and  E.  H.  Anderson,  "Detailed  Models  of  Piezoceramic  Actuation  of 
Beams,"  Journal  of  Intelligent  Material  Systems  and  Structures.  Vol.  1,  pp.  4-25,  Jan  1990. 


108 


39.  Obal,  M.  W.,  "Vibrations  Control  of  Flexible  Structures  using  Piezoelectric  Devices  as 
Sensors  and  Actuators,"  Doctorate  of  Philosophy  Dissertation  at  Georgia  Institute  of 
Technology,  Atlanta,  GA,  Sept  1986. 

40.  Hanagud,  S.,  Obal,  M.  W.,  and  M.  Meyyappa,  "Electronic  damping  techniques  and 
active  vibration  control,"  AIAA  Journal,  Vol.  23,  No.  5,  pp.  443-453,  May  1985. 

41.  Crawley,  E.  F.,and  J.  de  Luis,  "Use  of  Piezoelectric  Actuators  as  Elements  of 
Intelligent  Structures,"  AIAA  Journal.  Vol.  25,  No.  10,  pp.  1377-1385,  October  1987. 

42.  Hagood,  N.,  Crawley,  E.  F.,  Luis,  Javier  de,  and  Eric  H.  Anderson,  "Development  of 
Integrated  Components  for  Control  of  Intelligent  Structures,"  Proceedings  of  the  1988 
American  Control  Conference.  Vol.  3,  p  1890-6,  June  1988. 

43.  Wadin,  J.  R.,  Vemitron  Corporation  -  A  Guide  to  Modem  Piezoelectric  Ceramics. 
Sales  Manual,  Morgan  Matroc- Vemitron  Division,  El  Toro,  CA,  1994. 

44.  Hanagud,  S.,  Won,  C.  C.,  and  M.  W.  Obal,  "Optimal  Placement  of  Piezoceramic 
Sensors  and  Actuators,"  American  Control  Conference  3.  pp.  1884-1889,  1988. 

45.  Kulkarni,  G.  and  S.  V.  Hanagud,  "Modeling  Issues  in  the  Vibration  Control  with 
Piezoceramic  Actuators,"  Journal  of  Smart  Materials  and  Structures.  Vol  24,  pp.  7-17,  1991. 

46.  Shames,  I.  H.  and  C.  L.  Dym,  Energy  and  Finite  Element  Methods  in  Structural 
Mechanics.  Hemisphere  Publishing  Corporation,  Washington,  1985. 

47.  Wang,  C.,  Applied  Elasticity.  McGraw-Hill  Book  Company,  Inc.,  NY,  1953. 

48.  Buhariwala,  K.  J.,  "Dynamics  of  Viscoelastic  Structures,"  PhD  Disseration  University 
of  Toronto  -  Institute  For  Aerospace  Studies,  Toronto,  Canada,  1986. 

49.  Leissa,  A.  W.,  and  T.  H.  Young,  "Extensions  of  the  Ritz-Galerkin  Method  for  the 
Forced,  Damped  Vibrations  of  Structural  Elements,"  Proceedings  of  the  1984  Vibration 
Damping  Workshop,  pp.  EE-2  -  EE-20,  1984. 

50.  -  Flugge,  W.  (ed.),  Handbook  of  Engineering  Mechanics.  McGraw-Hill  Book  Company, 
New  York,  1962. 

51.  Lang,  D.  and  L.  Reithler,  "Active  Control  of  a  Composite  Plate:  From  the  Specific 
Problem  to  an  Optimized  Solution,"  Proceedings  of  1995  Smart  Structures  and  Materials 
Conference.  San  Diego,  CA,  Vol.  2443,  pp.  597-608,  Mar  1-2  1995. 

52.  Plass,  H.  J.,  Gaines,  J.,  H.,  and  C.  D.  Newson,  "Application  of  Reissner’s  Variational 
Principle  to  Cantilever  Plate  Deflection  and  Vibration  Problems,"  Journal  of  Applied 
Mechanics.  Vol  29,  #1,  pp.  127-135,  Mar  1962. 


109 


53.  Wolfram,  S.,  Mathematica  -  A  System  for  Doing  Mathematics  bv  Computer.  Addison- 
Wesley  Publishing  Co.,  Reading,  MA,  1988. 

54.  Bodden,  D.  S.,  and  J.  L.  Junkins,  "Eigenvalue  Optimization  Algorithms  for 
Structure/Control  Design  Iterations,"  AIAA  Journal  of  Guidance.  Control  and  Dynamics.  Vol  8 
No.  6  ,  pp  697-706,  Nov  -  Dec  1985. 

55.  Vanderplaats,  Miura  &  Associates,  Inc.,  Dot  Users  Manual  Version  4.00.  VMA 
Engineering,  1993. 

56.  Slater,  G.  L.,  "A  Disturbance  Model  for  the  Optimization  of 
Control/Structure  Interactions  for  Flexible  Dynamic  Systems,"  AIAA  Guidance,  Navigation 
and  Control  Conference,  Minneapolis,  NM,  August  15  -  17,  1988,  pp  57-63;  AIAA  Paper  88- 
4058-CP. 

57.  Valentino,  J.,  Fortran  for  Technologists  and  Engineers.  Holt,  Rinehart,  and  Winston, 
Chicago,  IL,  1985. 

58.  Carlsson,  L.  A.  and  R.  B.  Pipes,  Experimental  Characterization  of  Advanced 
Composite  Materials.  Prentice-Hall,  Inc.,  Englewood  Cliffs,  NJ,  1987. 


110 


A  ppendix  A 


The  information  in  this  appendix  will  show  the  development  of  how  the  kinetic,  strain, 
and  work  energy  expressions  are  integrated  by  parts.  In  the  energy  expressions  there  are  four 
different  integrations  to  be  performed  for  the  dependent  variables.  The  first  involves  only  one 
derivative  with  respect  to  x  on  the  variation: 


a 

f A(x ,  t)  5 B(x,  t)  ,  =  A(x,  t)  5 B{x,  t )  j 

J  o 

°  (158) 

-Ja(x,  t)  ,xbB{x,  t)  dx 
o 


This  integral  contributes  to  the  x  boundary  condition  and  the  equation  of  motion. 

The  second  integral  to  be  performed  involves  two  derivatives  with  respect  to  x  on  the 
variation. 
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This  integral  contributes  to  the  x  boundary  conditions  and  the  equation  of  motion. 

The  third  integral  to  be  performed  involves  one  derivative  with  respect  to  t  on  the 
variation. 


Ill 


(160) 


J fA(x*  t)  5 B{x,  t )  dxdt 

Cx  0 


Ja  ( x ,  t )  SB  ( x ,  t)  \  dx 


^2  a 

-  JJa  (x,  t)  5B  (x,  t)  dxdt 
u  o 


This  integral  contributes  to  both  the  initial  conditions  and  the  equation  of  motion. 

The  fourth  integral  to  be  performed  involves  one  derivative  with  respect  to  t  and  one 
with  respect  to  x  on  the  variation. 
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This  integral  contributes  to  both  the  initial  conditions  and  the  equation  of  motion. 
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A  ppendix  B 


The  optimization  program  is  divided  up  into  many  different  parts.  There  is  an  overall  shell 
program  which  calls  the  various  subroutines  to  perform  the  optimization.  Figure  31  shows  a 
flow  chart  diagraming  this  optimization  program. 


Figure  31  Optimization  Program  Flow  Chart 
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This  program  is: 


c  Optimization  Program 
c 

c  Define  Constants  Make  All  Variables  Double  Precision 
c 

DOUBLE  PRECISION  X(14),  OBJ,  XL(14),  XU(14),  G(2), 
*  WK(2000),  RPRM(200) 
c 

DIMENSION  IWK(2000),  IPRM(200) 
c 

c  Define  NRWK,  NRIWK 
c 

NRWK  =  2000 
NRIWK  =  2000 
c 

c  Zero  RPRM  And  IPRM 
c" 

DO  10  I  =  1,200 
RPRM(I)  =  0.0D0 
10  IPRM®  =  0 
c 

c  Define  Method,  NDV,  NCON 
c  Sequential  Quadratic  Programming  Method 
c 

METHOD  =  1 
c 

c  Design  Variables 
c 

NDV  =  14 
c 

c  2  Constraints 
c 

NCON  =  2 
c 

c  Define  Bounds  and  Initial  Design  For  Theta’s 
c 

DO  20  K  =  1,10 
X(K)  =  1.570796327D0 
c 

c  Lower  Bounds 
c 

XL(K)  =  -1.570796327D0 
c 
c 
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20  XU(K)  =  1.570796327D0 
c 

c  Define  Bounds  and  Initial  Design  For  Other  Variables 

c 

c 

X(1 1)  =  4  8. DO  - 

X(12)  =  16. DO 
X(  1 3)  =  .0000 1D0 
X(14)  =  .00001D0 
c 

c  Lower  Bounds 
c 

XL(ll)  =  1.25D0 
XL(12)  =  1.25D0 
XL(13)  =  0.0D0 
XL(14)  =  0.0D0 
c 

c  Upper  Bounds 
c 

XU(ll)  =  22.75D0 
XU(12)  =  22.75D0 
XU(13)  =  400000. DO 
XU(14)  =  400000. DO 
c 

c  Define  IPRINT,  MINMAX,  INFO 
c 

c  Print  Control 
c 

IPRINT  =  3 
c 

c  MINIMIZE 
c 

MINMAX  =  -1 
c 

c  Initialize  INFO  To  Zero 
c  - 

INFO  =  0 
c 

c  Optimize 
c 

100  CALL  DOT  (INFO,  METHOD,  IPRINT,  NDV,  NCON,  X,  XL,  XU,  OBJ, 
*  MINMAX,  G,  RPRM,  IPRM,  WK,  NRWK,  IWK,  NRIWK) 
c 

c  Finished 
c 

IF(INFO'.EQ.O)  STOP 
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c 


c  Evaluate  Objective  and  Constraint 
c 

Call  Subprog(OBJ,  X,  G) 
c 

c  Go  Continue  with  Optimization 
c 

GO  TO  100 
c 

end 

The  first  subroutine  called  is  the  optimization  program  DOT.  This  program  is  too  long  to  list 
in  this  document.  The  second  subroutine  call  evaluates  the  objective  function  and  the 
constraints.  This  program  is: 

SUBROUTINE  Subprog(OBJ,  X,  G) 
c 

c  Subroutine  to  Calculate  Objective  Function 
c 

c  Define  Constants  Make  All  Variables  Double  Precision 
c 

INTEGER  LDSF,LDMF,LDEVEC,N,K,NOUT,LDS,LDM,IROW(8),NA,NDIM 
c 

PARAMETER  (LDM=4,  LDMINV=4,  LDS=4,  LDT=4,  NCIM=4,  NCS=4,  NCT=4, 

*  NRIM=4,  NRS=4,  NRT=4,  N=4,  LDT1=4,  NCT1=2,  NRT1=4,  LDT2=4, 

*  NCT2=2,  NRT2=4,  LDT3=4,  NCT3=1,  NRT3=4,  LDT4=4,  NCT4=1,  NRT4=4, 

*  LDG=2,  NCG=2,  NRG=2,  LDC=2,  NCC=8,  NRC=2,  LDT5=2,  NCT5=8, 

*  NRT5=2,  LDB1=8,  NCB1=2,  NRB1=8,  LDT6=8,  NCT6=8,  NRT6=8,  LDSF=4, 

*  LDMF=4,  LDEVEC=4,  NCD=1,  NRD=1,  LDD=1,  NCH=1,  NRH=8,  LDH=8, 

*  NCB=8,  NRB=8,  LDB=8,  NCBD=8,  NRBD=8,  LDBD=8,  NCA=8,  NRA=8, 

*  LDA=8,  NRASYS=8,  NCASYS=8,  LDASYS=8,  NCW1=8,  NRW1=1,  LDW1=1, 

*  NCW2=8,  NRW2=1,  LDW2=1,  NCT8=8,  NRT8=1,  LDT8=1,  NCT11=8, 

*  NRT1 1=1,  LDT1 1=1,  NCT9=1,  NRT9=1,  LDT9=1,  NCXRMS=1,  NRXRMS=8, 

*  LDXRMS=8,  NRT10=1,  NCT10=1,  LDT10=1,  NRT12=2,  NCT12=8,  LDT12=2, 

-*  NRT13=2,  NCT13=8,  LDT13=2,  NRT14=2,  NCT14=2,  LDT14=2,  NRT15=2, 

*  NCT15=2,  LDT15=2,  NCW3=8,  NRW3=1,  LDW3=1) 
c 

DOUBLE  PRECISION  El,  E2,  nul2,  nu21,  gl2,  rho,  e31,  actloc, 

*  wd,  Clle,  Rf,  A  1(4),  A2(4),  A3(4),  A4(4),  dill,  dl  12,  dl  13,  dll, 

*  d  1 6 1 ,  dl62,  dl63,  dl64,  dl6,  d661,  d662,  d663,  d66, 

*  mfree(LDMF,N),  sfree(LDSF,N),  AMACH,  BETA(N),  mass(LDM,LDM), 

*  invmas(LDMINV,LDMINV),  stiff(LDS,NCS),  KD,  KS,  d31, 

*  GAIN(LDG,NCG),  xd,  Ls,  Asys(LDASYS,NCASYS),  B1(LDB1,NCB1), 

*  C(LDC,NCC),  H(8,l),  A(LDA,NCA),  TEMP(LDT,NCT),  TEMP  1  (LDT 1  ,NCT  1 ) , 

*  TEMP2(LDT2,NCT2),  TEMP3(LDT3,NCT3),  TEMP4(LDT4,NCT4), 
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*  TEMP5(LDT5,  NCT5),  TEMP6(LDT6,NCT6),  X(14),  G(2),  OBJ, 

*  BETAl(N),  B(LDB,NCB),  BD(LDBD,NCBD),  D(LDD,  NCD), 

*  TEMP7(LDH,  NCH),  W1(LDW1,NCW1),  W2(LDW2,NCW2),  TEMP8(LDT8,NCT8), 

*  ZETA,  TEMP9(LDT9,NCT9),  TEMP10(LDT10,  NCT10), 

*  TEMPI  1(LDT11,  NCT11),  TEMP12(LDT12,NCT12),  TEMP13(LDT13,  NCT13), 

*  thick, TEMP  1 4(LDT  14,  NCT14),  TEMP15(LDT15,  NCT15),  W3(LDW3,NCW3)  - 
c 

DOUBLE  COMPLEX  ALPHA(N),  EVAL(N),  EVEC(LDEVEC,N),  ALPHA  1(N), 

*  EVALl(N),  EVEC1(LDEVEC,N) 
c 

EXTERNAL  DMRRRR,  AMACH,  GPIGR,  DGVCRG,  UMACH,  DWRRRN 
c 

NA  =  8 
NDIM  =  8 
c 

c  Define  Composite  Constants 
c 

El  =  22691551. 44D0 
E2  =  1 429537. 875D0 
gl2  =  868816. 825D0 
nul2  =  .305616286D0 
nu21  =  .019253424D0 
rho  =  .052598707D0 
c 

c  Weighting  Matrices  for  the  Performance  Index 
c 

ZETA  =  .5D0 
c 

c  Location  of  Disturbance  Torque 
c 

xd  =  24. DO 
c 

c  Strength  of  the  Disturbance  Torque 
c 

-  D(  1,1)  =  1.D0 
c 

c  Define  PZT  Actuator  and  Sensor  Constants 
c 

wd  =  1.5D0 
Rf  =  l.D+7 
Cl  le  =  8.844D+06 
Ls  =  2.5D0 
actloc  =  .025D0 
d31  =  6.7323D-9 
e31  =  Cl le*d31 
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KS  =  -wd*actloc*d31*Clle*Rf 
KD  =  e31*wd*actloc 
c 

c  Evaluate  Composite  Constants 
c 

dill  =  2.54167D-6*El*Cos(X(l))**4.D0/(l.D0-nul2*nu21)  +  - 

*  2.54167D-6*El*Cos(X(10))**4.D0/(l. DO-nu  12*nu21)+ 

*  1.54167D-6*El*Cos(X(2))**4.D0/(l.D0-nul2*nu21)+ 

*  7.91667D-7*El*Cos(X(3))**4.D0/(l.D0-nul2*nu21)+ 

*  2.91667D-7*El*Cos(X(4))**4.D0/(l.D0-nul2*nu21)+ 

*  4. 1 6667D-8*El*Cos(X(5))**4.D0/(  1  .DO-nu  1 2*nu21)+ 

*  4.16667D-8*El*Cos(X(6))**4.D0/(l.D0-nul2*nu21  )+ 

*  2.91667D-7*El*Cos(X(7))**4.D0/(l. DO-nu  12*nu21)+ 

*  7.9 1 667D-7*El*Cos(X(8))**4.D0/(  1  .DO-nu  1 2*nu21  )+ 

*  1.54167D-6*El*Cos(X(9))**4.D0/(l.D0-nul2*nu21)+ 

*  0.00001 01 6666666666667D0*gl2*Cos(X(l))**2*Sin(X(l))**2.D0+ 

*  5.08333D-6*E2*nul2*Cos(X(l))**2.D0*Sin(X(l))**2.D0/(l.D0-nul2* 

*  nu21)+2.54167D-6*E2*Sin(X(l))**4.D0/(l.D0-nul2*nu21)+ 

*  0.00001 01 6666666666667D0*gl2*Cos(X(10))**2.D0*Sin(X(10))**2.D0+ 

*  5.08333D-6*E2*nul2*Cos(X(10))**2.D0* 

*  Sin(X(10))**2.D0/(  1. DO-nu  12*nu21)+ 

*  2.541 67D-6*E2*Sin(X(10))**4.D0/(l. DO-  nul2*nu21) 
c 

dl  12  =  6.16667D-6*gl2*Cos(X(2))**2.D0*Sin(X(2))**2.D0  + 

*  3.08333D-6*E2*nul2*Cos(X(2))**2.D0*Sin(X(2))**2.D0/(l.D0-nul2* 

*  nu2 1  )+l  .541 67D-6*E2*Sin(X(2))**4.D0/(  1  .DO-  nul2*nu21)  + 

*  3.16667D-6*gl2*Cos(X(3))**2.D0*Sin(X(3))**2.D0  + 

*  1.58333D-6*E2*nul2*Cos(X(3))**2.D0*Sin(X(3))**2.D0/(l.D0-nul2* 

*  nu21)+7.91667D-7*E2*Sin(X(3))**4.D0/(l.D0-  nul2*nu21)  + 

*  1.16667D-6*gl2*Cos(X(4))**2.D0*Sin(X(4))**2.D0  + 

*  5.83333D-7*E2*nul2*Cos(X(4))**2.D0*Sin(X(4))**2.D0/(l.D0-nul2* 

*  nu2 1  )+2. 9 1 667D-7*E2*Sin(X(4))* *4.D0/(  1  .DO-  nul2*nu21)  + 

*  1.66667D-7*gl2*Cos(X(5))**2.D0*Sin(X(5))**2.D0  + 

*  8.33333D-8*E2*nul2*Cos(X(5))**2.D0*Sin(X(5))**2.D0/(l.D0-nul2* 

-*  nu21)+4.16667D-8*E2*Sin(X(5))**4.D0/(l.D0-  nul2*nu21)  + 

*  1.66667D-7*gl2*Cos(X(6))**2.D0*Sin(X(6))**2.D0  + 

*  8.33333D-8*E2*nul2*Cos(X(6))**2.D0*Sin(X(6))**2.D0/(l.D0-nul2* 

*  nu21) 
c 

dl  13  =  4.16667D-8*E2*Sin(X(6))**4.D0/(l.D0-  nul2*nu21)  + 

*  1.16667D-6*gl2*Cos(X(7))**2.D0*Sin(X(7))**2.D0  + 

*  5.83333D-7*E2*nul2*Cos(X(7))**2.D0*Sin(X(7))**2.D0/(l.D0-nul2* 

*  nu21)+2.91667D-7*E2*Sin(X(7))**4.D0/(l.D0-  nul2*nu21)  + 

*  3.16667D-6*gl2*Cos(X(8))**2.D0*Sin(X(8))**2.D0  + 

*  1.58333D-6*E2*nul2*Cos(X(8))**2.D0*Sin(X(8))**2.D0/(l  .D0-nul2* 
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*  nu21)+7.91 667D-7*E2*Sin(X(8))**4.D0/(  1  .DO-  nul2*nu21)  + 

*  6.16667D-6*gl2*Cos(X(9))**2.D0*Sin(X(9))**2.D0  + 

*  3.08333D-6*E2*nul2*Cos(X(9))**2.D0*Sin(X(9))**2.D0/(l.D0-nul2* 

*  nu21  )+l  .541 67D-6*E2*Sin(X(9))**4.D0/(  1  .DO-  nu!2*nu21) 

dll  =  dill  +  dl  12  +  dl  13 

d  1 6 1  =  -5.08333D-6*gl  2*Cos(X(l))**3.D0*Sin(X(l ))  + 

*  2.54167D-6*E1  *Cos(X(l))**3.DO*Sin(X(l))/(l  .DO-  nul2*nu21)  - 

*  2.54167D-6*E2*nul2*Cos(X(l))**3.D0*Sin(X(l))/(l.D0-nul2*nu21)+ 

*  5.08333D-6*gl2*Cos(X(l))*Sin(X(l))**3.D0  - 

*  2.54167D-6*E2*Cos(X(l))*Sin(X(l))**3.D0/(l.D0-nul2*nu21)  + 

*  2.54167D-6*E2*nul2*Cos(X(l))*Sin(X(l))**3.D0/(l.D0-nul2*nu21)- 

*  5.08333D-6*gl2*Cos(X(10))**3.D0*Sin(X(10))  + 

*  2.54167D-6*El*Cos(X(10))**3.D0*Sin(X(10))/(l.D0-nul2*nu21)- 

*  2.541 67D-6*E2*nul 2*Cos(X(10))**3. DO*Sin(X(l 0))/(l.D0-nul 2*nu21)+ 

*  5.08333D-6*gl2*Cos(X(10))*Sin(X(10))**3.D0  - 

*  2.54167D-6*E2*Cos(X(10))*Sin(X(10))**3.D0/(l.D0-nul2*nu21)  + 

*  2.541 67D-6*E2*nul2*Cos(X(10))*Sin(X(  10))**3.D0/(1. D0-nul2*nu21)- 

*  3.08333D-6*gl2*Cos(X(2))**3.D0*Sin(X(2))  + 

*  1.54167D-6*El*Cos(X(2))**3.D0*Sin(X(2))/(l.D0-nul2*nu21)  - 

*  1.54167D-6*E2*nul2*Cos(X(2))**3.D0*Sin(X(2))/(l.D0-nul2*nu21)+ 

*  3.08333D-6*gl2*Cos(X(2))*Sin(X(2))**3.D0  - 

*  1.54167D-6*E2*Cos(X(2))*Sin(X(2))**3.D0/(l.D0-nul2*nu21) 

dl62  =  1 .54167D-6*E2*nul2*Cos(X(2))*Sin(X(2))**3.D0/(l.D0-nul2* 

*  nii21)-1.58333D-6*gl2*Cos(X(3))**3.D0*Sin(X(3))  + 

*  7.91667D-7*El*Cos(X(3))**3.D0*Sin(X(3))/(l.D0-  nul2*nu21)  - 

*  7.91667D-7*E2*nul2*Cos(X(3))**3.D0*Sin(X(3))/(l.D0-nul2*nu21)+ 

*  1.58333D-6*gl2*Cos(X(3))*Sin(X(3))**3.DO  - 

*  7.91667D-7*E2*Cos(X(3))*Sin(X(3))**3.D0/(l.D0-  nul2*nu21)  + 

*  7.91667D-7*E2*nul2*Cos(X(3))*Sin(X(3))**3.D0/(l  .D0-nul2*nu21)- 

*  5.83333D-7*gl2*Cos(X(4))**3.D0*Sin(X(4))  + 

*  2.91667D-7*El*Cos(X(4))**3.D0*Sin(X(4))/(l.D0-  nul2*nu21)  - 

-*  2.91667D-7*E2*nul2*Cos(X(4))**3.D0*Sin(X(4))/(l.D0-nul2*nu21)+ 

*  5.83333D-7*gl2*Cos(X(4))*Sin(X(4))**3.D0  - 

*  2.91667D-7*E2*Cos(X(4))*Sin(X(4))**3.D0/(l.D0-  nul2*nu21)  + 

*  2.91667D-7*E2*nul2*Cos(X(4))*Sin(X(4))**3.D0/(l.D0-nul2*nu21)- 

*  8.33333D-8*gl2*Cos(X(5))**3.DO*Sin(X(5))  + 

*  4.16667D-8*El*Cos(X(5))**3.D0*Sin(X(5))/(l.D0-  nul2*nu21)  - 

*  4.16667D-8*E2*nul2*Cos(X(5))**3.D0*Sin(X(5))/(l.D0-nul2*nu21)+ 

*  8.33333D-8*gl2*Cos(X(5))*Sin(X(5))**3.DO  - 

*  4.l6667D-8*E2*Cos(X(5))*Sin(X(5))**3.D0/(l.D0-  nul2*nu21)  + 

*  4.16667D-8*E2*nul2*Cos(X(5))*Sin(X(5))**3.D0/(l.D0-nul2*nu21) 
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d  1 63  =  -8.33333D-8*gl2*Cos(X(6))**3.D0*Sin(X(6))  + 

*  4.16667D-8*El*Cos(X(6))**3.D0*Sin(X(6))/(l.D0-  nul2*nu21)  - 

*  4.16667D-8*E2*nul2*Cos(X(6))**3.D0*Sin(X(6))/(l.D0-nul2*nu21)+ 

*  8.33333D-8*gl2*Cos(X(6))*Sin(X(6))**3.D0  - 

*  4.16667D-8*E2*Cos(X(6))*Sin(X(6))**3.D0/(l.D0-nul2*nu21)  + 

*  4.16667D-8*E2*nul2*Cos(X(6))*Sin(X(6))**3.D0/(l.D0-nul2*nu21)- 

*  5.83333D-7*gl2*Cos(X(7))**3.D0*Sin(X(7))  + 

*  2.91667D-7*El*Cos(X(7))**3.D0*Sin(X(7))/(l.D0-  nul2*nu21)  - 

*  2.91667D-7*E2*nul2*Cos(X(7))**3.D0*Sin(X(7))/(l.D0-nul2*nu21)+ 

*  5.83333D-7*gl2*Cos(X(7))*Sin(X(7))**3.D0  - 

*  2.91667D-7*E2*Cos(X(7))*Sin(X(7))**3.D0/(l.D0-  nul2*nu21)  + 

*  2.91667D-7*E2*nul2*Cos(X(7))*Sin(X(7))**3.D0/(l.D0-nul2*nu21)- 

*  1 ,58333D-6*gl  2*Cos(X(8))**3.D0*Sin(X(8))  + 

*  7.91667D-7*El*Cos(X(8))**3.D0*Sin(X(8))/(l.D0-  nul2*nu21)  - 

*  7.91667D-7*E2*nul2*Cos(X(8))**3.D0*Sin(X(8))/(l.D0-nul2*nu21)+ 

*  1.58333D-6*gl2*Cos(X(8))*Sin(X(8))**3.D0  - 

*  7.91667D-7*E2*Cos(X(8))*Sin(X(8))**3.D0/(l.D0-  nul2*nu21)  + 

*  7.91667D-7*E2*nul2*Cos(X(8))*Sin(X(8))**3.D0/(l.D0-nul2*nu21) 


dl64  =  -3.08333D-6*gl2*Cos(X(9))**3.D0*Sin(X(9))  + 

*  1.54167D-6*El*Cos(X(9))**3.D0*Sin(X(9))/(l.D0-  nul2*nu21)  - 

*  1.54167D-6*E2*nul2*Cos(X(9))**3.D0*Sin(X(9))/(l.D0-nul2*nu21)+ 

*  3.08333D-6*gl2*Cos(X(9))*Sin(X(9))**3.D0  - 

*  1.54167D-6*E2*Cos(X(9))*Sin(X(9))**3.D0/(l.D0-nul2*nu21)+ 

*  1.54167D-6*E2*nul2*Cos(X(9))*Sin(X(9))**3.D0/(l.D0-nul2*nu21) 

dl6  =  dl61  +  dl62  +  dl63  +  dl64 

d661  =  2.541 67D-6*gl2*Cos(X(l))**4.D0  + 

*  2.54167D-6*gl2*Cos(X(10))**4.D0  + 

*  1.54167D-6*gl2*Cos(X(2))**4.D0+7.91667D-7*gl2*Cos(X(3))**4.D0+ 

*  2.91667D-7*gl2*Cos(X(4))**4.D0+4.16667D-8*gl2*Cos(X(5))**4.D0+ 

*  4.16667D-8*gl2*Cos(X(6))**4.D0+2.91667D-7*gl2*Cos(X(7))**4.D0+ 

*  7.9 1667D-7*gl  2*Cos(X(8))**4.D0+ 1.541 67D-6*gl  2*Cos(X(9))**4.D0- 
-*  5.08333D-6*gl2*Cos(X(l))**2.D0*Sin(X(l))**2.D0  + 

*  2.54167D-6*El*Cos(X(l))**2.D0*Sin(X(l))**2.D0/(l.D0-nul2*nu21)+ 

*  2.54167D-6*E2*Cos(X(l))**2.D0*Sin(X(l))**2.D0/(l.D0-nul2*nu21)- 

*  5.08333D-6*E2*nul2*Cos(X(l))**2.D0* 

*  Sin(X(l))**2.D0/(l.D0-nul2*nu21)+2.54167D-6*gl2*Sin(X(l))**4.D0- 

*  5.08333D-6*gl2*Cos(X(10))**2.D0*Sin(X(10))**2.D0  + 

*  2.54167D-6*El*Cos(X(10))**2.D0*Sin(X(10))**2.D0/(l.D0-nul2*nu21)+ 

*  2.541 67D-6*E2*Cos(X(10))**2.D0*Sin(X(10))**2.D0/( l.D0-nul2*nu21  )- 

*  5.08333D-6*E2*nul2*Cos(X(10))**2.D0* 

*  Sin(X(10))**2.D0/(l.D0-nul2*nu21)+2.54167D-6*gl2* 

*  Sin(X(l0))*M.D0-3.08333D-6*gl2*Cos(X(2))**2.D0*Sin(X(2))**2.D0+ 
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*  1.54167D-6*El*Cos(X(2))**2.D0*Sin(X(2))**2.D0/(l.D0-  nul2*nu21) 

d662  =  1.54167D-6*E2*Cos(X(2))**2.D0* 

*  Sin(X(2))**2.D0/(l.D0-nul2*nu21)- 

*  3.08333D-6*E2*nul2*Cos(X(2))**2.D0*Sin(X(2))**2.D0/(l.D0-nul2* 

*  nu21)+1.54167D-6*gl2*Sin(X(2))**4.D0  - 

*  1.58333D-6*gl2*Cos(X(3))**2.D0*Sin(X(3))**2.D0  + 

*  7.91667D-7*El*Cos(X(3))**2.D0*Sin(X(3))**2.D0/(l. DO-nu  12*nu21)+ 

*  7.91 667D-7*E2*Cos(X(3))**2.D0*Sin(X(3))**2.D0/(l. DO-nu  12*nu21)- 

*  1.58333D-6*E2*nul2*Cos(X(3))**2.D0* 

*  Sin(X(3))**2.D0/(l.D0-nul2*nu21)+7.91667D-7*gl2*Sin(X(3))**4.D0- 

*  5.83333D-7*gl2*Cos(X(4))**2.D0*Sin(X(4))**2.D0  + 

*  2.9l667D-7*El*Cos(X(4))**2.D0*Sin(X(4))**2.D0/(l.D0-nul2*nu21)+ 

*  2.9l667D-7*E2*Cos(X(4))**2.D0*Sin(X(4))**2.D0/(l.D0-nul2*nu21)- 

*  5.83333D-7*E2*nul2*Cos(X(4))**2.D0* 

*  Sin(X(4))**2.D0/(l. DO-nu  12*nu21)+2.9l667D-7*gl2*Sin(X(4))**4.D0- 

*  8.33333D-8*gl2*Cos(X(5))**2.D0*Sin(X(5))**2.D0  + 

*  4.16667D-8*El*Cos(X(5))**2.D0*Sin(X(5))**2.D0/(l.D0-nul2*nu21)+ 

*  4.16667D-8*E2*Cos(X(5))**2.D0*Sin(X(5))**2.D0/(l.D0-nul2*nu21)- 

*  8.33333D-8*E2*nul2*Cos(X(5))**2.D0* 


*  Sin(X(5))**2.D0/(l. DO-nu  12*nu21)+4.16667D-8*gl2*Sin(X(5))**4.D0 
c 

d663  =  -8.33333D-8*gl2*Cos(X(6))**2.D0*Sin(X(6))**2.D0  + 

*  4.16667D-8*El*Cos(X(6))**2.D0*Sin(X(6))**2.D0/(l.D0-nul2*nu21)+ 

*  4.16667D-8*E2*Cos(X(6))**2.D0*Sin(X(6))**2.D0/(l.D0-nul2*nu21)- 

*  8.33333D-8*E2*nul2*Cos(X(6))**2.D0* 

*  Sin(X(6))**2.D0/(l. DO-nu  12*nu21)+4.16667D-8*gl2*Sin(X(6))**4.D0- 

*  5.83333D-7*gl2*Cos(X(7))**2.D0*Sin(X(7))**2.D0  + 

*  2.91667D-7*El*Cos(X(7))**2.D0*Sin(X(7))**2.D0/(l.D0-nul2*nu21)+ 

*  2.9 1 667D-7*E2*Cos(X(7))**2.D0*Sin(X(7))**2.D0/(  1  .DO-nu  1 2*nu21 )- 

*  5.83333D-7*E2*nul2*Cos(X(7))**2.D0* 

*  Sin(X(7))**2.D0/(l. DO-nu  12*nu21)+2.91667D-7*gl2*Sin(X(7))**4.D0- 

*  1.58333D-6*gl2*Cos(X(8))**2.D0*Sin(X(8))**2.D0  + 

*  7.9 1667D-7*E1  *Cos(X(8))**2.D0*Sin(X(8))**2.D0/(l. DO-nu  12*nu21)+ 
-*  7.9 1667D-7*E2*Cos(X(8))**2.D0*Sin(X(8))**2.D0/(l. DO-nu  12*nu21)- 

*  1.58333D-6*E2*nul2*Cos(X(8))**2.D0* 

*  Sin(X(8))**2.D0/(  1  .DO-nu  1 2*nu21  )+7.91667D-7*g  1 2*Sin(X(8))**4.D0- 

*  3.08333D-6*gl2*Cos(X(9))**2.D0*Sin(X(9))**2.D0  + 

*  1. 541 67D-6*El*Cos(X(9))**2.D0*Sin(X(9))**2.D0/(l. DO-nu  12*nu21)+ 

*  1.54167D-6*E2*Cos(X(9))**2.D0*Sin(X(9))**2.D0/(l.D0-nul2*nu21)- 

*  3.08333D-6*E2*nul2*Cos(X(9))**2.D0* 

*  Sin(X(9))**2.D0/(l, DO-nu  12*nu21)+l. 54 167D-6*gl2*Sin(X(9))**4.D0 
c 

d66  =  d661  +  d662  +  d663 
c 
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c 

c 

c 

c 

c 


Calculate  The  Mass  (mfree)  And  Stiffness  (sfree)  Matrices  For  The 
Free  Motion  System 

Mass  Matrix 

mfree(l,l)  =  1.08845D0*rho 
mfree(l,2)  =  0.D0 
mfree(l,3)  =  2.76282D0*rho 
mfree(l,4)  =  0.D0 
mfree(2,l)  =  mfree(l,2) 
mfree(2,2)  =  0.36287400**0 
mfree(2,3)  =  O.DO 
mfree(2,4)  =  0.921084D0*rho 
mfree(3,l)  =  mfree(l,3) 
mfree(3,2)  =  mfree(2,3) 
mfree(3,3)  =  7.20001D0*rho 
mfree(3,4)  =  O.DO 
mfree(4,l)  =  mfree(l,4) 
mfree(4,2)  =  mfree(2,4) 
mfree(4,3)  =  mfree(3,4) 
mfree(4,4)  =  2.40038D0*rho 


c 

c 

c 


c 

c 

c 

c 


Stiffness  Matrix 

sfree(U)  -  0.000880796902434194D0*dl  1 
sfree(l,2)  =  -0.00856736D0*dl6 
sfree(l,3)  =  0.001495287262745941DO*dl  1 
sfree(l,4)  =  -0.0228463D0*dl6 
sfree(2,l)  =  sfree(l,2) 

sfree(2,2)  =  0.0002935989674780654D0*dl  l+0.205617D0*d66 

sfree(2,3)  =  0.02284630648400314D0*dl6 

sfree(2,4)  =  0.00049842908758 19803D0*dl  l+0.349066D0*d66 

sfree(3,l)  =  sfree(l,3) 

sfree(3,2)  =  sfree(2,3) 

sfree(3,3)  =  0.0140927504389471  lDO*dl  1 

sfree(3,4)  =  O.DO 

sfree(4,l)  =  sfree(l,4) 

sfree(4,2)  =  sfree(2,4) 

sfree(4,3)  =  sfree(3,4) 

sfree(4,4)  =  0.004697583479649038D0*dl  1  +  0.822467D0*d66 

Call  a  Subroutine  Program  Which  Evaluates  The  Eigenvalues  and 
Eigenvectors  For  The  Free  Motion  System 

CALL  DGVCRG  (N,  sfree,  LDSF,  mfree,  LDMF,  ALPHA,  BETA,  EVEC, 
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*  LDEVEC) 

DO  10  K=1,N 

IF  (BETA(K)  .NE.  0.0)  THEN 
EVAL(K)  =  ALPHA(K)/BETA(K) 

ELSE 

EVAL(K)  =  AMACH(2) 

END  IF 

10  CONTINUE 

Al(l)  =  EVEC(  1 , 1  )/(EVEC(  1 , 1  )**2.D0+EVEC(2, 1  )**2.D0+ 

*  EVEC(3,1)**2.D0+EVEC(4,1)**2.D0)**.5D0 

A 1  (2)  =  EVEC(2, 1  )/(EVEC(  1 , 1  )**2.D0+EVEC(2,1  )**2.D0+ 

*  EVEC(3,1)**2.D0+EVEC(4,1)**2.D0)**.5D0 

A 1  (3)  =  EVEC(3, 1)/(EVEC(  1 , 1  )**2.D0+EVEC(2, 1  )**2.D0+ 
'  *  EVEC(3,1)**2.D0+EVEC(4,1)**2.D0)**.5D0 

A  1(4)  =  EVEC(4,1)/(EVEC(1,1)**2.D0+EVEC(2,1)**2.D0+ 

*  EVEC(3,1)**2.D0+EVEC(4,1)**2.D0)**.5D0 


A2(l)  =  EVEC(  1 ,2)/(EVEC(  1 ,2)**2.D0+EVEC(2,2)**2.D0+ 

*  EVEC(3,2)**2.D0+EVEC(4,2)**2.D0)**.5D0 

A2(2)  =  EVEC(2,2)/(EVEC(1,2)**2.D0+EVEC(2,2)**2.D0+ 

*  EVEC(3,2)**2.D0+EVEC(4,2)**2.D0)**.5D0 

A2(3)  =  EVEC(3,2)/(EVEC(1,2)**2.D0+EVEC(2,2)**2.D0+ 

*  EVEC(3,2)**2.D0+EVEC(4,2)**2.D0)**.5D0 

A2(4)  =  EVEC(4,2)/(EVEC(1,2)**2.D0+EVEC(2,2)**2.D0+ 

*  EVEC(3,2)**2.D0+EVEC(4,2)**2.D0)**.5D0 

-  A3(l)  =  EVEC(1,3)/(EVEC(1,3)**2.D0+EVEC(2,3)**2.D0+ 

*  EVEC(3,3)**2.D0+EVEC(4,3)**2.D0)**.5D0 

A3(2)  =  EVEC(2,3)/(EVEC(1,3)**2.D0+EVEC(2,3)**2.D0+ 

*  EVEC(3,3)**2.D0+EVEC(4,3)**2.D0)**.5D0 

A3(3)  =  EVEC(3,3)/(EVEC(1,3)**2.D0+EVEC(2,3)**2.D0+ 

*  EVEC(3,3)**2.D0+EVEC(4,3)**2.D0)**.5D0 

A3(4)  =  EVEC(4,3)/(EVEC(1,3)**2.D0+EVEC(2,3)**2.D0+ 

*  EVEC(3,3)**2.D0+EVEC(4,3)**2.D0)**.5D0 
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A4(l)  =  EVEC(1,4)/(EVEC(1,4)**2.D0+EVEC(2,4)**2.D0+ 

*  EVEC(3,4)**2.D0+EVEC(4,4)**2.D0)**.5D0 

A4(2)  =  EVEC(2,4)/(EVEC(1,4)**2.D0+EVEC(2,4)**2.D0+ 

*  EVEC(3,4)**2.D0+EVEC(4,4)**2.D0)**.5D0 

A4(3)  =  EVEC(3,4)/(EVEC(1,4)**2.D0+EVEC(2,4)**2.D0+ 

*  EVEC(3,4)**2.D0+EVEC(4,4)**2.D0)**.5D0 

A4(4)  =  EVEC(4,4)/(EVEC(1,4)**2.D0+EVEC(2,4)**2.D0+ 

*  EVEC(3,4)**2.D0+EVEC(4,4)**2.D0)**.5D0 


Calculate  The  Mass  And  Stiffness  Matrices  For  The  Forced  Motion 
System 

Stiffness  Matrix 


stiff(  1,1)  =  0.0008807969024341 94D0*Al(l)**2.D0*dll  + 

555  0.0002935989674780654D0*Al(2)**2.D0*dl  1  + 

*  0.00299057452549 1 88DO*A  1  (1  )*A1  (3)*d  1 1  + 

*  0.0140927504389471  lD0*Al(3)**2.D0*dl  1  + 

*  0.000996858 175 16396 1D0*A  1(2)*  A 1  (4)*d  1 1  + 

*  0.004697583479649034D0*Al(4)**2.D0*dl  1-0.01 71 347D0*A1(1)*A1(2) 

*  dl6+0.04569261296800626D0*Al(2)*Al(3)*dl6-0.0456926D0*Al(l)* 

*  Al(4)*dl6+0.205617D0*Al(2)**2.D0*d66  +  0.698132D0*A1(2)*A1(4)* 

*  d66  +0.822467D0*Al(4)**2.D0*d66 

stiff(l,2)  =  0.0008807969024341 95D0*  A 1  ( 1  )*A2(  1  )*d  1 1  + 

*  0.000293598967478065D0*Al(2)*A2(2)*dl  1  + 

*  0.00149528726274594D0*A2(l)*Al(3)*dl  1  + 

*  0.00149528726274594D0*Al(l)*A2(3)*dl  1  + 

*  0.0140927504389471  lD0*Al(3)*A2(3)*dl  1  + 

-*  0.00049842908758 19808D0*A2(2)*Al(4)*dl  1  + 

*  0.00049842908758 19808D0*Al(2)*A2(4)*dl  1  + 

*  0.004697583479649042D0*A1  (4)*A2(4)*dl  1  - 

*  0.0O856736D0*A2(l)*Al(2)*dl6-0.00856736D0*Al(l)*A2(2)*dl6  + 

*  0.022846306484003 1 3D0* A2(2)* A 1  (3)*dl 6  + 

*  0.022846306484003 13D0*Al(2)*A2(3)*dl6  - 

*  0.0228463D0*A2(l)*Al(4)*dl6  -  0.0228463D0*Al(l)*A2(4)*dl6  + 

*  0.2056 1 7D0*A  1  (2)*A2(2)*d66  +  0.349066D0*A2(2)*Al(4)*d66  + 

*  0.349066D0*Al(2)*A2(4)*d66  +  0.822467D0*Al(4)*A2(4)*d66 


stiff(l,3)  =  0.000880796902434195D0*Al(l)*A3(l)*dl  1  + 


*  0.000293598967478065D0*Al(2)*A3(2)*dll  + 

*  0.00149528726274594D0*A3(l)*Al(3)*dl  1  + 

*  0.00 1 495287 26274594D0* A  1(1)  *  A3(3)*d  1 1  + 

*  0.0140927504389471  lD0*Al(3)*A3(3)*dll  + 

*  0.00049842908758 19805D0*A3(2)*Al(4)*dl  1  + 

*  0.0004984290875819805D0*Al(2)*A3(4)*dll  +  - 

*  0.004697583479649042D0*Al(4)*A3(4)*dl  1  - 

*  0.00856736DO*A3(l)*Al(2)*dl6-0.00856736DO*Al(l)*A3(2)*dl6  + 

*  0.022846306484003 13D0*A3(2)*Al(3)*dl 6  + 

*  0.022846306484003 13D0*Al(2)*A3(3)*dl 6  - 

*  0.0228463D0*A3(l)*Al(4)*dl6  -  0.0228463D0*Al(l)*A3(4)*dl6  + 

*  0.2056 17D0*Al(2)*A3(2)*d66  +  0.349066D0*A3(2)*Al(4)*d66  + 

*  0.349066D0*Al(2)*A3(4)*d66  +  0.822467D0*Al(4)*A3(4)*d66 

stiff(l,4)  =  0.000880796902434 195D0*Al(l)*A4(l)*dl  1  + 

*  0.000293598967478065D0*Al(2)*A4(2)*dl  1  + 

*  0.00149528726274594D0*A4(l)*Al(3)*dl  1  + 

*  0.00 1 49528726274594D0*A  1  (1  )*A4(3)*dl  1  + 

*  0.0140927504389471  lD0*Al(3)*A4(3)*dl  1  + 

*  0.00049842908758 19805D0*A4(2)*Al(4)*dl  1  + 

*  0.00049842908758 19805D0*Al(2)*A4(4)*dl  1  + 

*  0.004697 583479649042D0*A  1  (4)*A4(4)*dl  1  - 

*  0.00856736D0*A4(l)*Al(2)*dl6-0.00856736D0*Al(l)*A4(2)*dl6  + 

*  0.022846306484003 13D0*A4(2)*Al(3)*dl  6  + 

*  0.022846306484003 13D0*Al(2)*A4(3)*dl  6  - 

*  0.O228463D0*A4(l)*Al(4)*dl6-0.0228463D0*Al(l)*A4(4)*dl6  + 

*  0.205617D0*Al(2)*A4(2)*d66  +  0.349066D0*A4(2)*Al(4)*d66  + 

*  0.349066D0*Al(2)*A4(4)*d66  +  0.822467D0*Al(4)*A4(4)*d66 

stiff(2,l)  =  stiff ( 1 ,2) 

stiff(2,2)  =  0.0008807969024341 94D0*A2(l)**2.D0*dll  + 

*  0.0002935989674780654D0*A2(2)**2.D0*dl  1  + 

*  0.00299057452549 188D0*A2(l)*A2(3)*dl  1  + 

-*  0.0140927504389471  lD0*A2(3)**2.D0*dl  1  + 

*  0. 000996858 175 163961  D0*A2(2)* A2(4)*d  1 1  + 

*  0.004697583479649034D0*A2(4)**2.D0*dl  1-0.01 71 347D0*A2(1)*A2(2)* 

*  d  1 6+0.0456926 1 296800626D0* A2(2)*A2(3)*d  1 6-0.0456926D0*  A2(  1  )* 

*  A2(4)*dl  6+0.2056 1 7D0* A2(2)**2.D0*d66  +  0.698 132D0*A2(2)*A2(4)* 

*  d66  +0.822467D0*A2(4)**2.D0*d66 

stiff(2,3)  =  0. 0008807969024341 95D0*A2(l)*A3(l)*dl  1  + 

*  0.000293598967478065D0*A2(2)*A3(2)*dl  1  + 

*  0.0014952S726274594D0*A3(l)*A2(3)*dl  1  + 

*  0.00149528726274594D0*A2(l)*A3(3)*dl  1  + 
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*  0.0140927504389471  lD0*A2(3)*A3(3)*dll  + 

*  0.00049842908758 19805D0*A3(2)*A2(4)*d  11  + 

*  0.00049842908758 1 9805D0*A2(2)*A3(4)*d  1 1  + 

*  0.004697583479649042D0*A2(4)*A3(4)*dl  1  - 

*  0.00856736D0*A3(l)*A2(2)*dl6  -  0.00856736D0*A2(l)*A3(2)*dl6  + 

*  0.022846306484003 13D0*A3(2)*A2(3)*d  16  +  -  - 

*  0.022846306484003 13D0*A2(2)*A3(3)*dl6  - 

*  0.0228463D0*A3(l)*A2(4)*dl6  -  0.0228463D0*A2(l)*A3(4)*dl6  + 

*  0.205617D0*A2(2)*A3(2)*d66  +  0.349066D0*A3(2)*A2(4)*d66  + 

*  0.349066D0*A2(2)*A3(4)*d66  +  0.822467D0*A2(4)*A3(4)*d66 

stiff(2,4)  =  0. 000880796902434 195D0*A2(l)*A4(l)*dl  1  + 

*  0.000293598967478065D0*A2(2)*A4(2)*dl  1  + 

*  0.00 1 49528726274594D0*A4(1  )* A2(3)*d  1 1  + 

*  0.00149528726274594D0*A2(l)*A4(3)*dl  1  + 

*  0.0140927504389471  lD0*A2(3)*A4(3)*dl  1  + 

*  0.00049842908758 19808D0*A4(2)*A2(4)*dl  1  + 

*  0.00049842908758 19808D0*A2(2)*A4(4)*dl  1  + 

*  0.004697583479649042DO*A2(4)*A4(4)*dl  1  - 

*  0.00856736D0*A4(l)*A2(2)*dl6  -  0.00856736D0*A2(l)*A4(2)*dl6  + 

*  0.0228463064840031 3D0*A4(2)*A2(3)*dl 6  + 

*  0.022846306484003 13D0*A2(2)*A4(3)*d  16  - 

*  0.0228463D0*A4(l)*A2(4)*dl6  -  0.0228463D0*A2(l)*A4(4)*dl6  + 

*  0.205617D0*A2(2)*A4(2)*d66  +  0.349066D0*A4(2)*A2(4)*d66  + 

*  0.349066D0*A2(2)*A4(4)*d66  +  0.822467D0*A2(4)*A4(4)*d66 


stiff(3,l)  =  stiff(l,3) 
stiff(3,2)  =  stiff(2,3) 

stiff(3,3)  =  0.000880796902434 194D0*A3(l)**2.D0*dl  1  + 

*  0.0002935989674780654D0*A3(2)**2.D0*dl  1  + 

*  0.00299057452549 188D0*A3(l)*A3(3)*dl  1  + 

*  0.0140927504389471  lD0*A3(3)**2.D0*dl  1  + 

-*  0.000996858 175 16396 lD0*A3(2)*A3(4)*dl  1  + 

*  0.004697583479649034D0*A3(4)**2.D0*dl  1-0.01 7 1347D0*A3(1)*A3(2)* 

*  dl  6+0. 0456926 1296800626D0*A3(2)*A3(3)*dl6-0.0456926D0*A3(l)* 

*  A3(4)*dl  6+0. 2056 17D0*A3(2)**2.D0*d66+0. 698 132D0*A3(2)*A3(4)*d66+ 

*  0.822467D0*A3(4)**2.D0*d66 

stiff(3,4)  =  0.0008807969024341 95D0*A3(  1  )*A4(  1  )*d  1 1  + 

*  0.000293598967478065D0*A3(2)*A4(2)*dl  1  + 

*  0.00149528726274594D0*A4(l)*A3(3)*dl  1  + 

*  0.00149528726274594D0*A3(l)*A4(3)*dl  1  + 

*  0.0140927504389471  lD0*A3(3)*A4(3)*dll  + 
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*  0.00049842908758 19808D0*A4(2)*A3(4)*dl  1  + 

*  0.00049842908758 19808D0*A3(2)*A4(4)*dl  1  + 

*  0.004697583479649042D0*A3(4)*A4(4)*dl  1  - 

*  0.00856736D0*A4(l)*A3(2)*dl6  -  0.00856736D0*A3(l)*A4(2)*dl6  + 

*  0.022846306484003 13D0*A4(2)*A3(3)*dl6  + 

*  0.022846306484003 1 3D0*  A3(2)*A4(3)*dl  6  -  • 

*  0.0228463D0*A4(l)*A3(4)*dl6  -  0.0228463D0*A3(l)*A4(4)*dl6  + 

*  0.2056 17D0*A3(2)*A4(2)*d66  +  0.349066D0*A4(2)*A3(4)*d66  + 

*  0.349066D0*A3(2)*A4(4)*d66  +  0.822467D0*A3(4)*A4(4)*d66 


stiff(4,l)  =  stiff(l,4) 


stiff(4,2)  =  stiff(2,4) 


stiff(4,3)  =  stiff(3,4) 

stiff(4,4)  =  0.000S80796902434194D0*A4(l)**2.D0*dl  1  + 

*  0.0002935989674780654D0*A4(2)**2.D0*dl  1  + 

*  0. 00299057452549 188D0*A4(l)*A4(3)*dll  + 

*  0.0140927504389471  lD0*A4(3)**2.D0*dl  1  + 

*  0.000996858 175 163961  DO*  A4(2)*A4(4)*dl  1  + 

*  0.004697583479649034D0*A4(4)**2.D0*dl  1-0.017 1347D0*A4(1)*A4(2)* 

*  dl6+0.04569261296800626D0*A4(2)*A4(3)*dl6-0.0456926D0*A4(l)* 

*  A4(4)*dl 6+0.2056 17D0* A4(2)**2.D0*d66+0.698 132D0* A4(2)*A4(4)*d66+ 

*  0. 822467* A4(4)**2.D0*d66 


Mass  Matrix 


mass(l,l)  =  1.08845D0*Al(l)**2.D0*rho+0.362874D0*Al(2)**2.D0*rho+ 

*  5.52564D0*Al(l)*Al(3)*rho  +  7.20001D0*Al(3)**2.D0*rho  + 

*  1.84217D0*Al(2)*Al(4)*rho  +  2.40038D0*Al(4)**2.D0*rho 

mass(l,2)  =  1 ,08845D0*A  1  ( 1  )*  A2(  1  )*rho+0.362874D0*A  1  (2)*A2(2)*rho+ 

*  2.76282D0*A2(l)*Al(3)*rho  +  2.76282D0*Al(l)*A2(3)*rho  + 

-*  7.20001D0*Al(3)*A2(3)*rho  +  0.921084D0*A2(2)*Al(4)*rho  + 

*  0.921084D0*Al(2)*A2(4)*rho  +  2.40038D0*Al(4)*A2(4)*rho 

mass(l,3)  =  1.08845D0*Al(l)*A3(l)*rho+0.362874D0*Al(2)*A3(2)*rho+ 

*  2.76282D0*A3(l)*Al(3)*rho  +  2.76282D0*Al(l)*A3(3)*rho  + 

*  7.20001D0*Al(3)*A3(3)*rho  +  0.921084D0*A3(2)*Al(4)*rho  + 

*  0.921084D0*Al(2)*A3(4)*rho  +  2.40038D0*Al(4)*A3(4)*rho 

mass(l,4)  =  1.08845D0*Al(l)*A4(l)*rho+0.362874D0*Al(2)*A4(2)*rho+ 

*  2.76282D0*A4(l)*Al(3)*rho  +  2.76282D0*Al(l)*A4(3)*rho  + 

*  7.2000 lD0*Al(3)*A4(3)*rho  +  0.921084D0*A4(2)*Al(4)*rho  + 


*  0.92 1 084D0*A1  (2)* A4(4)*rho  +  2.40038D0*Al(4)*A4(4)*rho 


mass(2,l)  =  mass(l,2) 

mass(2,2)  =  1.08845D0*A2(l)**2.D0*rho+0.362874D0*A2(2)**2.D0*rho+ 

*  5.52564D0*A2(l)*A2(3)*rho  +  7.20001D0*A2(3)**2.D0*rho  +  - 

*  1.84217D0*A2(2)*A2(4)*rho  +  2.40038D0*A2(4)**2.D0*rho 

mass(2,3)  =  1.08845D0*A2(l)*A3(l)*rho+0.362874D0*A2(2)*A3(2)*rho+ 

*  2.76282D0*A3(l)*A2(3)*rho  +  2.76282D0*A2(l)*A3(3)*rho  + 

*  7.2000  lD0*A2(3)*A3(3)*rho  +  0.921084D0*A3(2)*A2(4)*rho  + 

*  0.921 084D0*A2(2)*A3(4)*rho  +  2.40038D0*A2(4)*A3(4)*rho 

mass(2,4)  =  L08845D0*A2(l)*A4(l)*rho+0.362874D0*A2(2)*A4(2)*rho+ 

*  2.76282D0*A4(l)*A2(3)*rho  +  2.76282D0*A2(l)*A4(3)*rho  + 

*  7.20001  DO*  A2(3)*A4(3)*rho  +  0.921084D0*A4(2)*A2(4)*rho  + 

*  0.921084D0*A2(2)*A4(4)*rho  +  2.40038D0*A2(4)*A4(4)*rho 


mass(3,l)  =  mass(l,3) 


mass(3,2)  =  mass(2,3) 

mass(3,3)  =  1.08845D0*A3(l)**2.D0*rho+0.362874D0*A3(2)**2.D0*rho+ 

*  5.52564D0*A3(l)*A3(3)*rho  +  7.20001D0*A3(3)**2.D0*rho  + 

*  1.84217D0*A3(2)*A3(4)*rho  +  2.40038D0*A3(4)**2.D0*rho 


mass(3,4)  =  1.08845D0*A3(l)*A4(l)*rho+0.362874D0*A3(2)*A4(2)*rho+ 

*  2.76282D0*A4(l)*A3(3)*rho  +  2.76282D0*A3(l)*A4(3)*rho  + 

*  7.2000 lD0*A3(3)*A4(3)*rho  +  0.921084D0*A4(2)*A3(4)*rho  + 

*  0.921084D0*A3(2)*A4(4)*rho  +  2.40038D0*A3(4)*A4(4)*rho 


mass(4,l)  =  mass(l,4) 


mass(4,2)  =  mass(2,4) 


mass(4,3)  =  mass(3,4) 

mass  (4,4)  =  1.08845D0*A4(l)**2.D0*rho+0.362874D0*A4(2)**2.D0*rho+ 

*  5.52564D0*A4(l)*A4(3)*rho  +  7.20001D0*A4(3)**2.D0*rho  + 

*  1.84217D0*A4(2)*A4(4)*rho  +  2.40038D0*A4(4)**2.D0*rho 


Call  a  Subroutine  Program  Which  Evaluates  The  Eigenvalues  and 
Eigenvectors  For  The  Free  Motion  System 


CALL  DGVCRG  (N,  stiff,  LDS,  mass,  LDM,  ALPHA  1,  BETA1,  EVEC1, 


*  LDEVEC) 

DO  100  K=1,N 

IF  (BETAl(K)  .NE.  0.0)  THEN 
EVALl(K)  =  ALPHA  1(K)/B  ETA  1(K) 
ELSE 

EVALl(K)  =  AMACH(2) 

END  IF 
100  CONTINUE 


Formulate  The  Open  And  Closed  Loop  State-Space  Equations 
Formulate  The  Inverse  Of  the  Mass  Matrix 


CALL  DLINRG  (N,  mass,  LDM,  invmas,  LDMINV) 


Multiply  The  Inverse  Of  the  Mass  Matrix  By  The  Stiffness 

CALL  DMRRRR  (NRIM,  NCIM,  invmas,  LDMINV,  NRS,  NCS,  stiff,  LDS, 

*  NRT,  NCT,  TEMP,  LDT) 

TEMP1(1,1)  =  KD*(0.0654498D0*A  1  ( 1  )*Sin(0.0654498D0*X(  1 1 ))  + 

*  0.1309D0*Al(3)*Sin(0.1309D0*X(l  1))) 

TEMP  1(2,1)  =  KD*(0.0654498D0*A2(l)*Sin(0.0654498D0*X(l  1))  + 

*  0.l309D0*A2(3)*Sin(0.1309D0*X(l  1))) 

TEMP1(3,1)  =  KD*(0.0654498D0*A3(l)*Sin(0.0654498D0*X(ll))  + 

*  0. 1309D0*A3(3)*Sin(0.1309D0*X(l  1))) 

TEMP  1(4,1)  =  KD*(0.0654498D0*A4(l)*Sin(0.0654498D0*X(l  1))  + 

*  0.1309D0*A4(3)*Sin(0.1309D0*X(l  1))) 

TEMP  1(1, 2)  =  KD*(0.0654498D0*  A 1  (1  )*Sin(0.0654498D0*X(  1 2))  + 

■*  0. 1309D0*Al(3)*Sin(0. 1309D0*X(12))) 

TEMPI  (2,2)  =  KD*(0.0654498D0*A2(l)*Sin(0.0654498D0*X(12))  + 

*  0. 1309D0*A2(3)*Sin(0. 1309D0*X(12))) 

TEMPI  (3,2)  =  KD*(0.0654498D0*A3(l)*Sin(0.0654498D0*X(12))  + 

*  0.1309D0*A3(3)*Sin(0.1309D0*X(12))) 

TEMP  1(4,2)  =  KD*(0.0654498D0*A4(l)*Sin(0.0654498D0*X(12))  + 

*  0.1309D0*A4(3)*Sin(0. 1309D0*X(12))) 


CALL  DMRRRR  (NRIM,  NCIM,  invmas,  LDMINV,  NRT1,  NCT1,  TEMPI,  LDT1, 

*  NRT2,  NCT2,  TEMP2,  LDT2) 

TEMP3(1,1)  =  -0.5D0*A1(2)  +  0.5D0*Al(2)*Cos(0.0654498D0*xd)  - 

*  0.5D0*A1(4)  +  0.5D0*A1  (4)*Cos(0. 1 309D0*xd) 

TEMP3(2,1)  =  -0.5D0*A2(2)  +  0.5D0*A2(2)*Cos(0.0654498D0*xd)  - 

*  0.5D0*A2(4)  +  0.5D0*A2(4)*Cos(0. 1 309D0*xd) 

TEMP3(3,1)  =  -0.5D0*A3(2)  +  0.5D0*A3(2)*Cos(0.0654498D0*xd)  - 

*  0.5D0*A3(4)  +  0.5D0*A3(4)*Cos(0.1309D0*xd) 

TEMP3(4,1)  =  -0.5D0*A4(2)  +  0.5D0*A4(2)*Cos(0.0654498D0*xd)  - 

*  0.5D0*A4(4)  +  0.5D0*A4(4)*Cos(0. 1309D0*xd) 

CALL  DMRRRR  (NRIM,  NCIM,  invmas,  LDMINV,  NRT3,  NCT3,  TEMP3,  LDT3, 

*  NRT4,  NCT4,  TEMP4,  LDT4) 


Derive  The  State  Space  Matrices 

Asys  Matrix 

Asys(l,l)  =  O.DO 
Asys(l,2)  =  O.DO 
Asys(I,3)  =  O.DO 
Asys(l,4)  =  0.0 
Asys(l,5)  =  EDO 
Asys(l,6)  =  O.DO 
Asys(l,7)  =  O.DO 

Asys(l,8)  =  O.DO 
Asys(2,l)  =  O.DO 
Asys(2,2)  =  O.DO 
Asys(2,3)  =  O.DO 
Asys(2,4)  =  O.DO 
Asys(2,5)  =  O.DO 
Asys(2,6)  =  l.DO 
Asys(2,7)  =  O.DO 
Asys(2,8)  =  O.DO 
Asys(3,l)  =  O.DO 
Asys(3,2)  =  O.DO 
Asys(3,3)  =  O.DO 
Asys(3,4)  =  O.DO 
Asys(3,5)  =  O.DO 
Asys(3,6)  =  O.DO 
Asys(3,7)  =  l.DO 


Asys(3,8)  =  O.DO 
Asys(4,l)  =  O.DO 
Asys(4,2)  =  O.DO 
Asys(4,3)  =  O.DO 
Asys(4,4)  =  O.DO 
Asys(4,5)  =  O.DO 
Asys(4,6)  =  O.DO 
Asys(4,7)  =  O.DO 
Asys(4,8)  =  l.DO 
Asys(5,l)  =  -TEMP(U) 
Asys(5,2)  =  -TEMP(1,2) 
Asys(5,3)  =  -TEMP(1,3) 
Asys(5,4)  =  -TEMP(1,4) 
Asys(5,5)  =  O.DO 
Asys(5,6)  =  O.DO 
Asys(5,7)  =  O.DO 
Asys(5,8)  =  O.DO 
Asys(6,l)  =  -TEMP(2,1) 
Asys(6,2)  =  -TEMP(2,2) 
Asys(6,3)  =  -TEMP(2,3) 
Asys(6,4)  =  -TEMP(2,4) 
Asys(6,5)  =  O.DO 
Asys(6,6)  =  O.DO 
Asys(6,7)  =  O.DO 
Asys(6,8)  =  O.DO 
Asys(7,l)  =  -TEMP(3,1) 
Asys(7,2)  =  -TEMP(3,2) 
Asys(7,3)  =  -TEMP(3,3) 
Asys(7,4)  =  -TEMP(3,4) 
Asys(7,5)  =  O.DO 


Asys(7,6)  =  O.DO 
Asys(7,7)  =  O.DO 
Asys(7,8)  =  O.DO 
Asys(8,l)  =  -TEMP(4,1) 
Asys(8,2)  =  -TEMP(4,2) 
Asys(8,3)  =  -TEMP(4,3) 
Asys(8,4)  =  -TEMP(4,4) 
Asys(8,5)  =  O.DO 
Asys(8,6)  =  O.DO 
Asys(8,7)  =  O.DO 
Asys(8,8)  =  O.DO 


B  1(1,2)  =  O.DO 
B  1(2,1)  =  O.DO 
B  1(2,2)  =  O.DO 
B  1(3,1)  =  O.DO 
B  1(3,2)  =  O.DO 

Bl(4,l)  =  O.DO  -■ 

B  1(4,2)  =  O.DO 
B 1  (5,1)  =  TEMP2(1,1) 

B  1(5,2)  =  TEMP2(1,2) 

B  1(6,1)  =  TEMP2(2,1) 

B  1(6,2)  =  TEMP2(2,2) 

B  1(7,1)  =  TEMP2(3,1) 

B  1(7,2)  =  TEMP2(3,2) 

B  1(8,1)  =  TEMP2(4,1) 

B  1(8,2)  =  TEMP2(4,2) 
c 

c  Gain  Matrix 

c 

GAIN(1,1)  =  X(13) 

GAIN(1,2)  =  O.DO 
GAIN(2,1)  =  O.DO 
GAIN(2,2)  =  X(14) 
c 

c  C  Matrix 

c 

C(l,l)  =  O.DO 
C(l,2)  =  O.DO 
C(l,3)  =  O.DO 
C(l,4)  =  O.DO 

C(l,5)=  -,025*KS*(0.0654498D0*Al(l)*Sin(0.0654498D0*(X(l  1)  + 

*  Ls/2.D0))  +  0. 1 309D0*A  1  (3)*Sin(0. 1 309D0*(X(  1 1 )  +  Ls/2.D0))  - 

*  0.0654498D0*Al(l)*Sin(0.0654498D0*(X(l  1)  -  Ls/2.D0))  - 

*  0.1309D0*Al(3)*Sin(0.1309D0*(X(l  1)  -  Ls/2.D0))) 

C(l,6)=  -.025*KS*(0.0654498D0*A2(l)*Sin(0.0654498D0*(X(l  1)  + 
-*  LS/2.D0))  +  0.1 309D0*A2(3)*Sin(0.1309D0*(X(l  1)  +  Ls/2.D0))  - 

*  0.0654498D0*A2(l)*Sin(0.0654498D0*(X(l  1)  -  Ls/2.D0))  - 

*  0. 1 309D0* A2(3)*Sin(0. 1 309D0*(X(1 1 )  -  Ls/2.D0))) 

C(l,7)=  -.025*KS*(0.0654498D0*A3(l)*Sin(0.0654498D0*(X(l  1)  + 

*  LS/2.D0))  +  0.1309D0*A3(3)*Sin(0.1309D0:,:(X(l  1)  +  Ls/2.D0))  - 

*  0.0654498D0*A3(l)*Sin(0.0654498D0*(X(l  1)  -  Ls/2.D0))  - 

*  0.1309D0*A3(3)*Sin(0.1309D0*(X(l  1)  -  Ls/2.D0))) 

C(l,8)=  -.025*KS*(0.0654498D0*A4(l)*Sin(0.0654498D0*(X(ll)  + 

*  Ls/2.D0))  +  0.1309D0*A4(3)*Sin(0.1309D0*(X(l  1)  +  Ls/2.D0))  - 

*  0.0654498D0*A4(l)*Sin(0.0654498D0*(X(l  1)  -  Ls/2.D0))  - 

*  0. 1 309D0*A4(3)*Sin(0. 1 309D0*(X(  11)  -  Ls/2.D0))) 
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C(2,l)  =  O.DO 
C(2,2)  =  O.DO 
C(2,3)  =  O.DO 
C(2,4)  =  O.DO 

C(2,5)=  -,025*KS*(0.0654498D0*Al(l)*Sin(0.0654498D0*(X(12)  + 

*  Ls/2.D0))  +  0.1309D0*Al(3)*Sin(0.1309D0*(X(12)  +  Ls/2.D0))- 

*  0.0654498D0*Al(l)*Sin(0.0654498D0*(X(12)  -  Ls/2.D0))  - 

*  0.1309D0*Al(3)*Sin(0.1309D0*(X(12)  -  Ls/2.D0))) 

C(2,6)=  -.025*KS*(0.0654498D0*A2(l)*Sin(0.0654498D0*(X(12)  + 

*  Ls/2.D0))  +  0. 1 309D0*A2(3)*Sin(0. 1 309D0*(X(  1 2)  +  Ls/2.D0))  - 

*  0.0654498D0*A2(  1  )*Sin(0.0654498D0*(X(  1 2)  -  Ls/2.D0))  - 

*  0.1309D0*A2(3)*Sin(0.1309D0*(X(12)  -  Ls/2.D0))) 

C(2,7)=  -.025*KS*(0.0654498D0*A3(l)*Sin(0.0654498D0*(X(12)  + 

*  Ls/2.D0))  +  0.1309D0*A3(3)*Sin(0.1309D0*(X(12)  +  Ls/2.D0))  - 

*  0.0654498D0*A3(l)*Sin(0.0654498D0*(X(12)  -  Ls/2.D0))  - 

*  0.1309D0*A3(3)*Sin(0.1309D0*(X(12)  -  Ls/2.D0))) 

C(2,8)=  -.025*KS*(0.0654498D0*A4(l)*Sin(0.0654498D0*(X(12)  + 

'  *  Ls/2.D0))  +  0.1309D0*A4(3)*Sin(0.1309D0*(X(12)  +  Ls/2.D0))  - 

*  0.0654498D0*A4(l)*Sin(0.0654498D0*(X(12)  -  Ls/2.D0))  - 

*  0. 1 309D0*A4(3)*Sin(0. 1 309D0*(X(1 2)  -  Ls/2.D0))) 
c 

c  H  Matrix 


H(l,l)  =  O.DO 
H(2,l)  =  O.DO 
H(3,l)  =  O.DO 
H(4,l)  =  O.DO 
H(5,l)  =  TEMP4(1,1) 
H(6,l)  =  TEMP4(2,1) 
H(7,l)  =  TEMP4(3,1) 
H(8,l)  =  TEMP4(4,1) 


c  Closed  Loop  Matrix 
c 

-  CALL  DMRRRR  (NRG,  NCG,  GAIN,  LDG,  NRC,  NCC,  C,  LDC,  NRT5,  NCT5, 

*  TEMP5,  LDT5) 
c 

CALL  DMRRRR  (NRB1,  NCB1,  Bl,  LDB1,  NRT5,  NCT5,  TEMP5,  LDT5, 

*  NRT6,  NCT6,  TEMP6,  LDT6) 
c 

A(  1,1)  =  Asys(l,l)  -  TEMP6(1,1) 

A(l,2)  =  Asys(l,2)  -  TEMP6(1,2) 

A(l,3)  =  Asys(l,3)  -  TEMP6(1,3) 

A(l,4)  =  Asys(l,4)  -  TEMP6(1,4) 

A(l,5)  =  Asys(l,5)  -  TEMP6(1,5) 
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A(l,6)  =  Asys(l,6)  -  TEMP6(1,6) 
A(l,7)  =  Asys(l,7)  -  TEMP6(1,7) 
A(l,8)  =  Asys(l,8)  -  TEMP6(1,8) 
A(2,l)  =  Asys(2,l)  -  TEMP6(2,1) 
A(2,2)  =  Asys(2,2)  -  TEMP6(2,2) 
A(2,3)  =  Asys(2,3)  -  TEMP6(2,3) 
A(2,4)  =  Asys(2,4)  -  TEMP6(2,4) 
A(2,5)  =  Asys(2,5)  -  TEMP6(2,5) 
A(2,6)  =  Asys(2,6)  -  TEMP6(2,6) 
A(2,7)  =  Asys(2,7)  -  TEMP6(2,7) 
A(2,8)  =  Asys(2,8)  -  TEMP6(2,8) 
A(3,l)  =  Asys(3,l)  -  TEMP6(3,1) 
A(3,2)  =  Asys(3,2)  -  TEMP6(3,2) 
A(3,3)  =  Asys(3,3)  -  TEMP6(3,3) 
A(3,4)  =  Asys(3,4)  -  TEMP6(3,4) 
A(3,5)  -  Asys(3,5)  -  TEMP6(3,5) 
A(3,6)  =  Asys(3,6)  -  TEMP6(3,6) 
A(3,7)  =  Asys(3,7)  -  TEMP6(3,7) 
A(3,8)  =  Asys(3,8)  -  TEMP6(3,8) 
A(4,l)  =  Asys(4,l)  -  TEMP6(4,1) 
A(4,2)  =  Asys(4,2)  -  TEMP6(4,2) 
A(4,3)  =  Asys(4,3)  -  TEMP6(4,3) 
A(4,4)  =  Asys(4,4)  -  TEMP6(4,4) 
A(4,5)  =  Asys(4,5)  -  TEMP6(4,5) 
A(4,6)  =  Asys(4,6)  -  TEMP6(4,6) 
A(4,7)  =  Asys(4,7)  -  TEMP6(4,7) 
A(4,8)  =  Asys(4,8)  -  TEMP6(4,8) 
A(5,l)  =  Asys(5,l)  -  TEMP6(5,1) 
A(5,2)  =  Asys(5,2)  -  TEMP6(5,2) 
A(5,3)  =  Asys(5,3)  -  TEMP6(5,3) 
A(5,4)  =  Asys(5,4)  -  TEMP6(5,4) 
A(5,5)  =  Asys(5,5)  -  TEMP6(5,5) 
A(5,6)  =  Asys(5,6)  -  TEMP6(5,6) 


A(5,7)  =  Asys(5,7)  -  TEMP6(5,7) 
A(5,8)  =  Asys(5,8)  -  TEMP6(5,8) 
A(6,l)  =  Asys(6,l)  -  TEMP6(6,1) 
A(6,2)  =  Asys(6,2)  -  TEMP6(6,2) 
A(6,3)  =  Asys(6,3)  -  TEMP6(6,3) 
A(6,4)  =  Asys(6,4)  -  TEMP6(6,4) 
A(6,5)  =  Asys(6,5)  -  TEMP6(6,5) 
A(6,6)  =  Asys(6,6)  -  TEMP6(6,6) 
A(6,7)  =  Asys(6,7)  -  TEMP6(6,7) 
A(6,8)  =  Asys(6,8)  -  TEMP6(6,8) 
A(7,l)  =  Asys(7,l)  -  TEMP6(7,1) 
Ml, 2)  =  Asys(7,2)  -  TEMP6(7,2) 
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A(7,3)  =  Asys(7,3)  -  TEMP6(7,3) 

A(7,4)  =  Asys(7,4)  -  TEMP6(7,4) 

A(7,5)  =  Asys(7,5)  -  TEMP6(7,5) 

A(7,6)  -  Asys(7,6)  -  TEMP6(7,6) 

A(7,7)  =  Asys(7,7)  -  TEMP6(7,7) 

A(7,8)  =  Asys(7,8)  -  TEMP6(7,8) 

A(8,l)  =  Asys(8,l)  -  TEMP6(8,1) 

A(8,2)  =  Asys(8,2)  -  TEMP6(8,2) 

A(8,3)  =  Asys(8,3)  -  TEMP6(8,3) 

A(8,4)  =  Asys(8,4)  -  TEMP6(8,4) 

A(8,5)  =  Asys(8,5)  -  TEMP6(8,5) 

A(8,6)  =  Asys(8,6)  -  TEMP6(8,6) 

A(8,7)  =  Asys(8,7)  -  TEMP6(8,7) 

A(8,8)  =  Asys(8,8)  -  TEMP6(8,8) 
c 

c  Form  the  Constant  Term  in  the  Lyap  Equation 
c 

CALL  DMRRRR  (NRH,  NCH,  H,  LDH,  NRD,  NCD,  D,  LDD,  NRH,  NCH, 

*  TEMP7,  LDH) 
c 

CALL  DMXYTF  (NRH,  NCH,  TEMP7,  LDH,  NRH,  NCH,  H,  LDH,  NRB, 

*  NCB,  B,  LDB) 

c 

c  Calculate  the  Lyap  Equation 

c 

Call  lyap(A,  B,  NA,  IROW,  NDIM) 
c 

c  Calculate  the  Performance  Index 

c 

Wl(l,l)  -  Al(l)  +  A  1(2)  +  2.D0*A1(3)  +  2.D0*A1(4) 

Wl(l,2)  =  A2(l)  +  A2(2)  +  2.D0*A2(3)  +  2.D0*A2(4) 

Wl(l,3)  =  A3(l)  +  A3(2)  +  2.D0*A3(3)  +  2.D0*A3(4) 

Wl(l,4)  =  A4(l)  +  A4(2)  +  2.D0*A4(3)  +  2.D0*A4(4) 

Wl(l,5)  -  0.D0 

-  Wl(l,6)  =  0.D0 
Wl(l,7)  =  0.D0 
Wl(l,8)  =  0.D0 
c 

W2(l,l)  =  Al(l)  +  2.D0*A1(3) 

W2(l,2)  =  A2(l)  +  2.D0*A2(3) 

W2(l,3)  =  A3(l)  +  2.D0*A3(3) 

W2(l,4)  =  A4(l)  +  2.D0*A4(3) 

W2(l,5)  =  0.D0 
W2(l,6)  =  0.D0 
W2(l,7)  =  0.D0 
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W2(l,8)  =  O.DO 
c 

W3(l,l)  =  Wl(l,l)  -  W2(l,l) 

W3(l,2)  =  W 1  ( 1 ,2)  -  W2(l,2) 

W3(l,3)  =  Wl(l,3)  -  W2(l,3) 

W3(l,4)  =  Wl(l,4)  -  W2(l,4)  ■■ 

W3(l,5)  =  Wl(l,5)  -  W2(l,5) 

W3(l,6)  =  Wl(l,6)  -  W2(l,6) 

W3(l,7)  =  Wl(l,7)  -  W2(l ,7) 

W3(l,8)  =  Wl(l,8)  -  W2(l,8) 
c 

CALL  DMRRRR  (NRW3,  NCW3,  W3,  LDW3,  NRB,  NCB,  B,  LDB, 

*  NRT8,  NCT8,  TEMP8,  LDT8) 
c 

CALL  DMXYTF  (NRT8,  NCT8,  TEMP8,  LDT8,  NRW3,  NCW3,  W3,  LDW3, 

*  NRT9,  NCT9,  TEMP9,  LDT9) 
c 

CALL  DMRRRR  (NRW2,  NCW2,  W2,  LDW2,  NRB,  NCB,  B,  LDB, 

*  NRT11,  NCT11,  TEMPI  1,  LDT11) 
c 

CALL  DMXYTF  (NRT11,  NCT11,  TEMP11,  LDT11,  NRW2,  NCW2,  W2, 

*  LDW2,  NRT10,  NCT10,  TEMPIO,  LDTIO) 
c 

OBJ  =  ZETA*TEMP9(  1,1)  +  (1-ZETA)*TEMP  10(1,1) 
c 

c  CALCULATE  THE  CONSTRAINTS 
c 

CALL  DMRRRR  (NRG,  NCG,  GAIN,  LDG,  NRC,  NCC,  C,  LDC, 

*  NRT12,  NCT12,  TEMP  12,  LDT12) 
c 

CALL  DMRRRR  (NRT12,  NCT12,  TEMP  12,  LDT12,  NRB,  NCB,  B,  LDB, 

*  NRT13,  NCT13,  TEMP  13,  LDT13) 
c 

CALL  DMXYTF  (NRT13,  NCT13,  TEMP  13,  LDT13,  NRC,  NCC,  C,  LDC, 

*  NRT14,  NCT14,  TEMP  14,  LDT14) 

c 

CALL  DMXYTF  (NRT14,  NCT14,  TEMP  14,  LDT14,  NRG,  NCG,  GAIN, 

*  LDG,  NRT15,  NCT15,  TEMP  15,  LDT15) 
c 

G(l)  =  TEMP  15(1,1)/1 60000. DO  -  1.D0 
c 

G(2)  =  TEMPI 5(2, 2)/l 60000.D0  -  1.D0 
c 

end 
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This  program  which  calculates  the  objective  function  and  constraints  also  calls  a  subroutine 
which  calculates  the  Lyapunov  of  the  function.  This  subroutine  is: 


SUBROUTINE  LYAP(A,B,NA,IROW,NDIM) 
implicit  double  precision(a-h,o-z) 

DOUBLE  PRECISION  A(NDIM,NDIM),B(NDIM,NDIM) 
DOUBLE  PRECISION  AMAT(1200,1200),RHS(1200),X(20,20) 
DOUBLE  PRECISION  WKAREA(  10000) 

INTEGER  IROW(NDIM) 

IA=1200 

DO  1  I=1,NA*NA 
DO  1  J=1,NA*NA 
AMAT(IJ)=0.0D0 
1  CONTINUE 
NUN=0 
DO  10  1=1,  NA 
DO  10  J=I,NA 
NUN=NUN+1 
RHS(NUN)=-B(I,J) 

NEQ=0 

DO  20  K=1,NA 
DO  20  L=K,NA 
NEQ=NEQ+1 
DO  30  M=1,NA 
DO  30  N=1,NA 
IF(K.EQ.M)THEN 

IF(I.EQ.N.AND.J.EQ.L)THEN 
AMAT(NEQ,NUN)=AMAT(NEQ,NUN)+A(M,N) 
END  IF 

IF(I.EQ.L.AND.J.EQ.N.AND.I.NEJ)THEN 
AMAT(NEQ,NUN)=AMAT(NEQ,NUN)+A(M,N) 
END  IF 
END  IF 

IF(L.EQ.M)THEN 

IF(I.EQ.N.AND.J.EQ.K)THEN 
AMAT(NEQ,NUN)=AMAT(NEQ,NUN)+A(M,N) 
END  IF 

IF(I.EQ.K.AND.J.EQ.N.AND.I.NE.J)THEN 
AMAT(NEQ,NUN)=AMAT(NEQ,NUN)+A(M,N) 
END  IF 
END  IF 

30  CONTINUE 
20  CONTINUE 
10  CONTINUE 

CALL  LUDCMP(AMAT,NEQ,IA,IROW,D) 

CALL  LUB  KSB  ( AM AT,NEQ,IA,IRO W,RHS) 
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c 


CALL  LEQT 1  F( AMAT,  1  ,NEQ,IA,RHS,0,WKAREA,IER) 
NEQ=0 

DO  31  1=1, NA 
DO  31  J=I,NA 
NEQ=NEQ+1 
B  (I,  J)=RHS  (NEQ) 

B(J,I)=B(I,J) 

31  CONTINUE 
RETURN 
END 

subroutine  lubksb(a,n,np,indx,b) 
implicit  double  precision(a-h,o-z) 
double  precision  a(np,np),b(n) 
integer  indx(np) 
ii=0 

do  12  i=l,n 
ll=indx(i) 
sum=b(ll) 
b(ll)=b(i) 
if  (ii.ne.O)then 
do  1 1  j=ii,i-l 
sum=sum-a(i,j)*bQ 

1 1  continue 

else  if  (sum.ne.O.)  then 
ii=i 
endif 
b(i)=sum 

12  continue 

do  14  i=n,  1,-1 
sum=b(i) 
if(i.lt.n)then 
do  13  j=i+l,n 
sum=sum-a(i,j)*b(j) 

13  continue 
endif 

b(i)=sum/a(i,i) 

14  continue 
return 
end 

subroutine  ludcmp(a,n,np,indx,d) 
implicit  double  precision(a-h,o-z) 
parameter  (nmax=500,tiny=l  ,0d-20) 
double  precision  a(np,np),vv(nmax) 
integer  indx(np) 
d=l. 

do  12  i=l,n 


138 


aamax=0.d0 
do  1 1  j=l,n 

if  (abs(a(i,j)).gt.aamax)  aamax=abs(a(i,j)) 

1 1  continue 

if  (aamax.eq.O.dO)  call  exit(l) 
vv(i)=l./aamax 

12  continue 
do  19  j=l,n 

if  (j.gt.l)  then 
do  14  i=l,j-l 
sum=a(i,j) 
if  (i.gt.l)then 
do  13  k=l,i-l 
sum=sum-a(i,k)*a(k,j) 

13  continue 
a(i,j)=sum 

endif 

14  continue 
endif 
aamax=0. 
do  16  i=j,n 

sum=a(i,j) 
if  (j.gt.l)then 
do  15  k=l,j-l 
sum=sum-a(i,k)*a(k,j) 

15  continue 
a(i,j)=sum 

endif 

dum=vv(i)*abs(sum) 
if  (dum.ge.aamax)  then 
imax=i 
aamax=dum 
endif 

16  continue 

if  (j.ne.imax)then 
do  17  k=l,n 
dum=a(imax,k) 
a(imax,k)=a(j,k) 
a(j,k)=dum 

17  continue 
d=-d 

vv(imax)=vv(j) 

endif 

indx(j)=imax 

if(j.ne.n)then 

if(a(jd).eq.0.d0  )a(j,j)=tiny 
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dum=l./a(j,j) 
do  18  i=j+l,n 
a(i,j)=a(i,j)*dum 

18  continue 
endif 

19  continue 

if(a(n,n).eq.0.d0)  a(n,n)=tiny 

return 

end 

subroutine  exit(iflag) 
if  (iflag  .ne.  0)  then 

print  *,  'singular  matrix.' 
stop 
endif 
return 
end 
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Appendix  C 


This  appendix  contains  the  raw  optimization  data  generated  from  the  optimization  program  for  each  starting 
location.  •  -  - 


Disturbance  Force  1 

Starting  Actuator  Location  8/16 

Starting  Gain  Multiplier  .00001 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

1.57079 

1.03-1.56 

ply  angle 

1.48352 

.519-1.47 

ply  angle 

1.39626 

0.88-1.39 

loc  1/loc  2 

8/16 

8,  19.2 

loc  1/loc  2 

8/16 

7.45,  22.7 

loc  1/loc  2 

’  8/16 

8,  22.7 

gain  1/2 

.00001 

0,  .0695 

gain  1/2 

.00001 

.29,  .12 

gain  1/2 

.00001 

.56,  .11 

objective 

4853.87 

objective 

467 

objective 

226 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

1.30899 

.78-1.3 

ply  angle 

1.22173 

.687-1.2 

ply  angle 

1.13446 

.286-1.12 

loc  1/loc  2 

8/16 

7.57,  22.7 

loc  1/loc  2 

8/16 

6.43,  22.7 

loc  1/loc  2 

8/16 

6.88,  22.7 

gain  1/2 

.00001 

.413,  1.2 

gain  1/2 

.00001 

.34,  .124 

gain  1/2 

.00001 

.929,  .127 

objective 

191 

objective 

171 

objective 

124 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

1.04719 

.42-1.04 

ply  angle 

0.95993 

.496-. 953 

ply  angle 

0.87266 

.51 1-.867 

loc  1/loc  2 

8/16 

7.29,  22.7 

loc  1/loc  2 

8/16 

7.52,  22.7 

loc  1/loc  2 

8/16 

7.7,  22.7 

gain  1/2 

.00001 

1.1,  .12 

gain  1/2 

.00001 

1.25,  .12 

gain  1/2 

.00001 

1.43,  .11 

objective 

138.4 

objective 

152 

objective 

158.8 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

0.78539 

.54-. 78 

ply  angle 

0.69813 

.592-. 696 

ply  angle 

0.61086 

.476-. 61 

loc  1/loc  2 

8/16 

7.6,  22.7 

loc  1/loc  2 

8/16 

6.7,  22.7 

loc  1/loc  2 

8/16 

4.9,  22.7 

gain  1/2 

.00001 

1.33,  1.14 

gain  1/2 

.00001 

.933,  .129 

gain  1/2 

.00001 

.743,  .14 

objective 

149.7 

objective 

130.5 

objective 

112.3 

Finish 

Start 

Finish 

Start 

Finish 

141 


ply  angle  0.52359 


loc  1/loc  2 


gain  1/2 


objective 


8/16 


.00001 


Start 


ply  angle  0.26179 


loc  1/loc  2 


gain  1/2 


objective 


8/16 


.00001 


.423-. 522 


5.03,  22.7 


.644,  .137 


115 


Finish 


.263-. 366 


5.24,  2.27 


.66.  .137 


106.3 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


0.43633  ,365-.435  ply  angle 


8/16  4.9,  22.7  loc  1/loc  2 


.00001  .669,  .138  gain  1/2 


109.3  objective 


Start  Finish 


0.17453  .181-. 59  ply  angle 


8/16  4.98,  22.7  loc  1/loc  2 


.00001  .50,  .134  gain  1/2 


121.9  objective 


0.34906  .34-. 35 


8/16  5.07,  22.7 


.00001  .676,  .138 


■  106.7 


Start  Finish 


0.08726  .097-. 298 


6,  22.7 


.00001  1.06, .14 


93.4 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle  - 


loc  1/loc  2 


gain  1/2 


objective 


Start 


1.57079 


8/16 


.0001 


Start 


1.30899 


8/16 


.0001 


Start 


1.04719 


8/16 


.0001 


Finish 


.383-1.54 


8,  22.7 


.529,  .113 


192 


Finish 


.56-1.3 


5,  22.7 


.55,  .135 


138 


Finish 


.67-1.04 


5.76,  22.7 


1.1,  .14 


159 


Finish 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


Start 


1.48352 


8/16 


.0001 


Start 


1.22173 


8/16 


.0001 


Start 


0.95993 


8/16 


.0001 


Finish 


1.18-1.48 


8,  20.6 


.0001,  .87 


4664 


Finish 


.304-1.2 


4.84,  22.7 


.73,  .14 


118.4 


Finish 


.644-.955 


5.77,  22.7 


1.1,  .14 


141 


Finish 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


Start 


1.39626 


8/16 


.0001 


Start 


1.13446 


8/16 


.0001 


Start 


0.87266 


8/16 


.0001 


Finish 


.878-1.4 


5.48,  22.7 


.78,  .14 


486 


Finish 


.691-1.13 


5.76,  22.7 


1.1,  .14 


183.5 


Finish 


.611-.87 


5.77,  22.7 


1.1,  .14 


125 


Finish 


ply  angle 

0.78539 

.586-/78 

ply  angle 

0.69813 

.24-. 69 

ply  angle 

0.61086 

.28-. 61 

loc  1/loc  2 

8/16 

5.78,  22.7 

loc  1/loc  2 

8/16 

6.01,  22.3 

loc  1/loc  2 

8/16 

5.98,  22.3 

gain  1/2 

.0001 

1.1,  .14 

gain  1/2 

.0001 

1.2,  .134 

gain  1/2 

.0001 

1.17,  .135 

objective 


ply  angle  1.30899  .286-1.3  ply  angle  1.22173  .31-1.21  ply  angle  1.13446  .45-1.12 


loc  1/loc  2 


sain  1/2 


5.9,  22.7 


143 


loc  1/loc  2 


aain  1/2 


6.45,  22.7 


ply  angle 


loc  1/loc  2 


sain  1/2 


0.52359  .226-. 52  ply  angle  0.43633  .325-.436  ply  angle  0.34906  .347-.36S 


6.3,  22.7 


ply  angle  1.57079  1.24-1.43  ply  angle  1.48352  .195-1.48  ply  angle  1.39626  .28-1.39 


loc  1/loc  2 


aain  1/2 


5.29,  22.7 


.852,  .14 


ply  angle  0.26179  .259-.347  ply  angle  0.17453  .18-.344  ply  angle  0.08726  .0935-.47 


loc  1/loc  2 


.sain  1/2 


6.1,  22.7 


ply  angle  1.57079  .866-1.55  ply  angle  1.48352  .353-1.48  ply  angle  1.39626  .349-1.39 


loc  1/loc  2 


sain  1/2 


8.12,  17.5 


145 


sain  1/2 


ply  angle  1.57079  1.38-1.4  ply  angle  1.48352  .378-1.48  ply  angle  1.39626  .346-1.39 


loc  1/loc  2 


sain  1/2 


8/16  6.09,  22.7 


ply  angle  1.30899  .321-1.29  ply  angle  1.22173  .34-1.16  ply  angle  1.13446  .352-1 


loc  1/loc  2 


sain  1/2 


8/16  6.53,  22.7 


8/16  7.05,  22.7  loc  1/loc  2 


8/16  7.15,  22.7 


1  5.83,  .143 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


1.04719 


8/16 


1 


Start 


0.78539 


8/16 


Start 


0.52359 


8/16 


1 


Start 


0.26179 


8/16 


1 


.35-. 74 


7.24,  22.7 


8.24,  .141 


61.9 


Finish 


.349-. 591 


7.19,  22.7 


6.57,  .143 


61.9 


Finish 


.338-.426 


7.22,  22.7 


7.36,  .142 


61.9 


Finish 


.32-.3S3 


7.24,  22.7 


7.93,  .142 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


0.95993 


8/16 


Start 


0.69813 


8/16 


Start 


0.43633 


8/16 


1 


Start 


0.17453 


8/16 


.352-. 669 


7.24,  22.7 


8.02,  .142 


61.9 


Finish 


.352-.494 


7.22,  22.7 


7.58,  .142 


61.76 


Finish 


,436-.444 


7.77,  22.7 


.8,  0 


2142 


Finish 


.267-.356 


7.21,  22.7 


7.03,  .142 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


ply  angle 


loc  1/loc  2 


gain  1/2 


0.87266 


8/16 


1 


Start 


0.61086 


8/16 


1 


Start 


0.34906 


‘8/16 


1 


Start 


0.08726 


8/16 


1 


.351-.565 


7.24,  22.7 


8.07,  .142 


61.8 


Finish 


.351-. 45 


7.22,  22.7 


7.47,  .142 


61.9 


Finish 


,339-.352 


7.22,  22.7 


7.43,  .142 


61.8 


Finish 


-1.57-.35 


7.13,  22.7 


5.28,  .142 


63.4 


Start 


ply  angle  .  1.57079 


loc  1/loc  2 


gain  1/2 


objective 


Starting  Gain  Multiplier  10 


Finish 


1.4 


8/16  7.39,  22.7 


10  8.86,  .12 


26875 


Finish 


ply  angle 


loc  1/loc  2 


gain  1/2 


objective 


Start  Finish 


1.48352  .344-1.48  ply  angle 


8/16  7.08,  22.7  loc  1/loc  2 


10  4.8,  .143  gain  1/2 


81.9  objective 


Finish 


Start  Finish 


1.39626  .372-1.39 


8/16  7.24,  22.7 


10  8.81,  .142 


83.3 


Start  Finish 


ply  angle 

1.30899 

.349-1.32 

ply  angle 

1.22173 

.333-1.21 

ply  angle 

1.13446 

.35-1.1 

loc  1/loc  2 

8/16 

7.25,  22.7 

loc  1/loc  2 

8/16 

7.24,  22.7 

loc  1/loc  2 

8/16 

7.24,  22.7 

gain  1/2 

10 

8.46,  .142 

gain  1/2 

10 

8.47,  .142 

gain  1/2 

10 

8.33,  .142 

objective 

66.2 

objective 

66.45 

objective 

66 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

1.04719 

.34-1.0 

ply  angle 

0.95993 

.33-. 93 

ply  angle 

0.87266 

.35-. 85 

loc  1/loc  2 

8/16 

7.24,  22.7 

loc  1  /loc  2 

8/16 

7.24,  22.7 

loc  1/loc  2 

8/16 

7.23,  22.7 

gain  1/2 

10 

8.23,  .142 

gain  1/2 

10 

8.2,  .142 

gain  1/2 

10 

8.17,  .142 

objective 

65 

objective 

64.2 

objective 

63.6 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

0.78539 

.351-. 761 

ply  angle 

0.69813 

.354-.679 

ply  angle 

0.61086 

.31-. 6 

loc  1/loc  2 

8/16 

7.24,  22.7 

loc  1/loc  2 

8/16 

7.24,  22.7 

loc  1/loc  2 

8/16 

7.25,  22.7 

gain  1/2 

10 

8.18,  .142 

gain  1/2 

10 

8.21,  .142 

gain  1/2 

10 

8.4,  .142 

objective 

63.1 

objective 

62.7 

objective 

62.45 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

0.52359 

.364-.517 

ply  angle 

0.43633 

.33-. 435 

ply  angle 

0.34906 

.348-. 36 

loc  1/loc  2 

8/16 

7.23,  22.7 

loc  1/loc  2 

8/16 

7.25,  22.7 

loc  1/loc  2 

8/16 

7.23,  22.7 

gain  1/2 

10 

7.95,  .142 

gain  1/2 

10 

8.46,  .142 

gain  1/2 

10 

7.93,  .142 

objective 

62.2 

objective 

61.8 

objective 

61.5 

Start 

Finish 

Start 

Finish 

Start 

Finish 

ply  angle 

0.26179 

.263-. 33 

ply  angle 

0.17453 

.183-.387 

ply  angle 

0.08726 

.288-. 356 

loc  1/loc  2 

8/16. 

7.46,  22.7 

loc  1/loc  2 

8/16 

7.21,  22.7 

loc  1/loc  2 

8/16 

7.24,  22.7 

gain  1/2 

10 

7.67,  0 

gain  1/2 

10 

7.12,  .142 

gain  1/2 

10 

8.09,  .142 

objective 

3059 

objective 

62.3 

objective 

61.67 
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Vila 


Capt  Douglas  DeHart  he 

graduated  from  James  B.  Conant  High  School  and  then  went  to  the  University  of  Wisconsin  at 
Madison  on  a  four  -  year  ROTC  Scholarship.  In  1985.  he  graduated  from  the  University  of 
Wisconsin  at  Madison  with  a  B.S.  in  Engineering  Mechanics  and  was  commissioned  a  2nd  Lt 
in  the  Air  Force.  He  then  received  an  education  delay  from  AFIT/CI  to  obtain  a  Master’s 
Degree.  In  Dec.  1986,  Capt  DeHart  graduated  from  the  University  of  Wisconsin  at  Madison 
with  a  M.S.  in  Engineering  Mechanics  emphasizing  in  Aerospace  Engineering.  He  entered 
active  duty  in  1987  at  Edwards  AFB,  California  and  was  assigned  to  the  Rocket  Propulsion 
Laboratory  as  a  Space  Systems  Technology  Project  Analyst.  During  his  stay  at  Edwards  AFB, 
the  laboratory  changed  its  name  and  expanded  its  mission.  It  became  the  Air  Force 
Astronautics  Laboratory  and  became  the  Air  Force's  main  laboratory  in  space  structure 
research.  At  the  Astronautics  Laboratory,  Capt  DeHart  managed  both  contracts  and  conducted 
research  in  the  area  of  smart  structures  or  embedded  sensors  arid  actuators.  Based  upon  his 
research  at  the  Astronautics  Laboratory,  he  applied  for  and  was  accepted  for  the  PhD  program 
at  the  Air  Force  Institute  of  Technology  in  April  1990.  He  started  his  studies  in  October  1994 
and-in  April  1994,  he  left  AF1T  and  was  assigned  to  the  Air  Force  Phillips  Laboratory,  one  of 
the  new  Air  Force  "super"  laboratories.  While  finishing  up  his  research.  Capt  DeHart  also 
manage  contracts  in  the  area  of  multi-functional  structures  which  is  the  integration  of  the 
satellite  housekeeping  functions  into  the  load  bearing  structural  shear  panels  which  make  up 
the  spacecraft  bus.  Currently,  he  is  the  branch  chief  of  the  Applied  Composite  Branch.  He  is 
in  charge  of  12  engineers  and  3  technicians.  They  conduct  research  and  manage  contracts  to 
promote  the  use  of  composite  materials  in  both  spacecraft  and  launch  vehicles. 
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